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1.0 INTRODUCTION 

This report represents an extension and refinement of programs and 
techniques developed for determining the vertical distributions of various 
atmospheric constituents (aerosols, ozone, and neutral atmospheric density 
i.e. Rayleigh scattering air molecules) from radiance measurements of 
the earth's horizon (Newell and Gray, 1972). Surveyed in the previously 
reported study were the techniques required for determining selected 
constituent distributions and the observability conditions and instrument 
characteristics necessary for the successful measurement inversion. This 
report extends the radiance inversion technique to the determination of 
aerospl physical characteristic distributions. In the previous study 
stratospheric aerosols were found to exhibit gross extinction features so 
that they could be observed by the proposed technique. Because of increasing 
interest in these particles as products of various pollution sources, it was 
decided to investigate further (Grayet al, 19 72) the potential of the horizon 
inversion technique for determining other physical characteristics of the 
stratospheric and mesospheric aerosols besides total extinction. Thus an 
analysis of the capability of the technique to deduce information about the 
particulate size distribution and index of refraction was undertaken. 

These investigations were centered around the refinement and extended 
development of the techniques previously employed. Here the state vector 
in the inversion computer code was extended to include the desired aerosol 
characteristics. The radiative transfer simulation was refined to include 
arbitrary (noncoplanar) azimuth angles and opaque (thick) clouds. These 
refinements were tested for accuracy against other less efficient codes. 
To supplement these activities an efficient empirical aerosol model was 
developed from data generated by a more complex Mie code computation 
of the aerosol optical properties. With the empirical aerosol model, 
quantities such as the partial derivatives of the angular scattering function 
with respect to the real part of the index of refraction can be easily computed. 
These refined techniques and aerosol models were then used to simulate 
an aerosol inversion procedure for obtaining the vertical distributions of 
aerosols and their gross physical properties. 
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The measurement technique used here is to scan the earth's sunlit 
horizon as shown in Figure 1.0-1 with a multiwavelength photometer and 
to record the radiant intensities measured at a predetermined sequence of 
angular positions. These measured intensities are then processed as shown 
in the flow chart Figure 1.0-2. _ .... _ •. 

The estimation procedure begins with an input of the geometry which 
includes primarily the sun direction and the measurement platform altitude. 
The geometric data along with an a priori estimate of the atmospheric 
state (constituent densities and aerosol characteristics) is used to compute 
a theoretical prediction of the intensity profile and its partial derivatives 
with respect to state elements. This is done with a radiative transfer 
simulation for the first selected wavelength and scan angle. The theoretical 
partial derivatives, theoretical intensities, and the measured intensity are 
then fed to the filter equations where an optimal estimate of the atmospheric 
state is produced based upon the difference between measured and predicted 
intensities, the state covariance, and the measurement noise. The optimal 
state estimate is then fed back through the intensity calculations and a 
new intensity is computed and compared to the next measurement. The 
loop continues until all wavelength channels and scan angles have been 
sampled. The result is an optimal estimate, in a minimum variance sense, 
of the atmospheric state including aerosol number density, size distribution 
parameters, and index of refraction. 

Section 2 of the report discusses the aerosol inversion capability, of 
the scattered sunlight experiment, stellar occultation inversion, and radiance 
effects from a horizon profile parameter study. The aerosol inversion 
simulations show the relative invertability of the various aerosol 
characteristics under a range of conditions including different altitudes 
and assumed a priori uncertainties. The efficacies of several wavelength 
bands are compared. Three aerosol size ranges are considered, (0.01-0.1*), 
<0.1 -1.0*) and (1.0-10*). The effects of changes in variables such as season, 
latitude, satellite altitude, solar zenith and azimuth angles, and albedo on 
horizon profiles are investigated. 
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Section 3 details the filter and radiative transfer model development 
and inversion procedure refinement work that was carried out to make 
aerosol inversion possible. Included are discussions of the augmented state 
filter, the refined radiative transfer simulation, and the empirical aerosol 
models. 

Section 4 provides a review of our conclusions and recommendations 
on the invertibility of aerosols, model development and satellite data 
requirements. 

Section 5 appends supplementary information. 

Appendix 1. A parameter study of horizon profile intensity variability 
as a function of season, latitude, albedo, sun azimuths, and zenith and 
wavelength is illustrated. 

Appendix 2. The regression formulation of the aerosol model used 
for the size ranges (0*01-0. 1*/), (0.1-1. 0a/) and (1.0-10.0/u). 

Appendix 3. The multiple scattering code C.K.W. is discussed. 
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Fig. 1.0-2 Atmospheric Data Inversion Procedure 
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2.0 RADIANCE INVERSION AND SIMULATION 

2.1 AEROSOLS - IMPLICATIONS OF INVERSION 

2.1.1 Introduction 

Using a Kalman-Bucy filter to invert horizon profile data produces 
both an estimate of the atmospheric state and a covariance matrix describing 
the accuracy of the estimates. By propagating the covariance matrix through 
a proposed measurement schedule it is possible to predetermine which 
parameters will or will not be affected by the measurements, and what 
linear combinations of. parameters will be well known. A function called 
sensitivity has been defined which is a measure of the amount of information 
obtained about a particular state element in a measurement. The sensitivity 
curves are computed very rapidly, but fail to predict quantitatively the 
final variances and covariances of the parameters which are obtained by a 
covariance propagation. They do, however, qualitatively predict which 
parameters will have small variances at the end of a measurement schedule, 
and so they are used extensively in the following discussion. 

A sensitivity analysis was applied to the problem of determining which 

aerosol parameters are observable from a horizon scan experiment. It 

has been demonstrated (New;ell and Gray, 1972) that aerosol extinction per 

unit volume is readily recovered from limb scan data at all altitudes between 

-3 

10 and 85 km with an idealized instrument having a noise value of 10 

times the maximum horizon signal and between 10 and 100 km for a noise 

-4 

value of 10 . The sensitivity study, therefore, is aimed at determining 

the observability of four different parameters, namely aerosol number 

density (p), size distribution parameter (a), real index of refraction (n), 

and imaginary (complex) index of refraction (n'). r The size distribution 

parameter (a) is the radius exponent where the power law size distribution 

-ct 

is of the form n(r) = Ar . The results are broken into three groups 
representing three assumed aerosol size ranges: 0.01-0.1/u, O.l-l.O/y, and 
1.0-10 v. For these size ranges solar zenith angles were varied from +60° 
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to -30° with resulting scattering angles varying from approximately 50° 
to 140°. Constituent results based on these variable conditions are developed 
in Figures 2.1-1 to 2.1-12. Analysis shows how the invertability of the 
parameters is affected by 1) the geometry 2) the instrument noise, and 3) 
the statistics of the initial estimates of the aerosol parameters. This 
analysis was conducted for spherical aerosols and nonspherical aerosol 
shape factors which have a spherical equivalent. 

2.1.2 (0.1 -1.0//) Aerosol Inversion Simulations 

For this aerosol size range (0.1 -1.0//), as is also the case for the 

other size ranges, the aerosol number density is initially adjusted so that 

extinction at 5500 A is equal to the extinction used as the standard state 

(Newell and Gray, 1972) based on Elterman, 1966andl970, given the initial 

estimate that a = 3, n = 1.5, and n' = 0. Because of the proportionality 

between extinction and cross section for a given number density there will 

be an increased number density for the (0.01-0.1//) range over the (0.1 -1.0//) 

range and conversely a decrease for the (1.0-10//) range. For the range 

3 3 

(0.1 -1.0//) thenumber density runs from, approximately 10 /cm at ground 
to 10 ^ / cm^ at 100 km. 

For all of the three aerosol size ranges the instrument white noise 
"7 2 

was assumed to be 10 /u watts/cm (RMS). It has been shown previously 
(Newell and Gray, 1972) that varying the noise simply varies the maximum 
altitude from which information can be retrieved. The maximum useful 
altitude can be approximated by observing the intensity versus altitude 
curves of Section 2.3 and noting at what altitude the intensity reaches the 
noise level i.e. S/N = 1. At a few kilometers above this point the initial 
estimates of constituents are as good as the filtered values and no significant 
changes are affected in the. estimates. 

Sensitivity curves were produced for a number of different scattering 
angles and wavelengths with an instrument having a finite field of view. 
Sensitivity is a function derived from the filter equations, which indicates 
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the relative information content from constituent to constituent and altitude 
to altitude for every measurement condition i.e. (wavelength and tangent 
height). Given a measurement, the information content for a particular 
constituent and altitude (state element) is reflected in a decrease in the 
variance associated with that element and it is natural to look at the 
covariance update equation for the definition of sensitivity: 

P (m+1) = P(m) - K (m+1) B (m+1) P (m) (2.1-1) 

the term of interest here is K(m+l)B(m+l)P(m), which is the decrease in 
the covariance for a particular measurement condition. Normalizing by 
P(m) yields the sensitivity matrix K(m+l)B(m+l). ,Of particular interest , 
in the sensitivity matrix are the diagonal elements which indicate the relative 
decrease in variance of the elements of the state vector. To isolate the 
effects of the measurement vector on the sensitivity, the covariance is 
assumed to maintain its initial value. With this assumption the sensitivity 
of the jth element of the state vector is 



In choosing an instrument's wavelengths it is desirable to have as 
many linearly independent measurements as possible. That is, each 
measurement should be sensitive to a different combination of parameters. 
The functions which relate extinction and phase function to aerosol 
parameters change slowly with respect to wavelength. For this reason a 
wide spread of wavelengths over the region 2500 to 7000 A is best suited 
for the inversion (see Figure 2.1-1 or 2.1-4). Thus the wavelengths 3000, 
4000, 5000, 6000, and 7000 A, which span the visible region with a minimum 
of wavelength channels necessary for the inversion, have been arbitrarily 
selected as working values. In this case the actual wavelengths chosen 
are not as important as the separation between the wavelengths and in an 
actual experiment more wavelengths would be run for greater accuracy. " 
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The total scattered sunlight radiance at the horizon is a composite 
of the radiance contributed by aerosols and Rayleigh scatterers. For each 
scattering constituent the radiance contribution is proportional to the product 
of its scattering cross section and angular scattering function (phase 
function). If the scattering cross .sections for both aerosol and Rayleigh 
scatterers are equal then the relative contribution of each constituent to 
the horizon radiance depends only upon the relative values of the phase 
functions at a given scattering angle. For most scattering angles away 
from the forward region the aerosol phase function will be smaller than 
the Rayleigh function thus reducing the relative contribution of aerosols to 
the radiance, however this effect can be partly offset by an increase in the 
relative scattering cross section of aerosols at longer wavelengths. 

Both' the effects of changes in wavelength on cross section and changes 
in scattering angle on the phase function can.be observed by examining the 
sensitivity curves of Figures 2.1-2 and 2.1-4. In Figure 2.1-2 it can be 
seen that the aerosol number density sensitivity ’’AEROSOL" increases at 
longer wavelengths while the "NEUTRAL" molecule scattering decreases. 
This is the cross-section effect.. Figure 2.1-4 differs only by a change in 
solar zenith from Figure 2.1-2. The scattering angle for Figure 2.1-4 is 
less favorable (i.e. less forward, «110° versus w 50°) for aerosols. 
Therefore, the aerosol sensitivity is lowered relative to Figure 2.1-2 while 
the neutral molecule scattering sensitivity increases. 

The most serious problem encountered is in the uncertainty of the 
initial statistics associated with the initial estimate. The updating of 
estimates by the filter inversion routine depends in part on the RMS values 
of the initial estimate. If a parameter is well known it will not receive 
much of an update and its variance will decrease relatively slowly. This 
is particularly critical in the case of the aerosol parameters where there 
is a large sensitivity to aerosol extinction per unit volume, however, the 
particular aerosol parameters which receive these updates are determined 
partly by their initial variances. 
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Sensitivity computations in the (0.1 -1.0//) size range were made with 
initial RMS error values of 600% for (p), 1 for (a), 0.05 for (n), and 0.005 
for (n 1 ) along with realistic estimates of neutral density and ozone (Malchow, 
1971) RMS values as a function of altitude (Figure 2.1-3). Curves run 
with these values showed a large sensitivity to aerosol number density 
and only a small additional sensitivity to (n 1 ) below 30 km. Other runs 
with 100% RMS error on (p) (Figure 2.1-1) showed all parameters except 
(n) to be readily observable at all altitudes. Thus varying the initial RMS 
estimates varies the final RMS estimates and the sensitivity, and until 
realistic data are collected about these parameters the covariance matrices 
associated with the parameter estimates will be unreliable. 

The sensitivity curves indicate a potential problem with regard to ob- 
servability, Since the wavelength dependence of sensitivity is basically the 
same for both number density(p) "AEROSOL" and size distribution a) "ALPHA 11 
there is the possibility that no two measurements are linearly independent in 
these two parameters, therefore making it impossible to estimate either in 
one , This similarity in wavelength dependence occurs because increasing 
either (a) or (p) has the effect of lowering the wavelength dependence of 
the total received signal by increasing the effective aerosol to Rayleigh 
scattering ratio. A complete propagation of the covariance matrix was 
made to determine if this was the case. The results indeed showed a 
high correlation (0.6 to 0.8) between the two parameters but also showed a 
significant factor of five decrease in variance for both parameters in a 
coarse inversion (only three wavelengths and a four kilometer separation 
of tangent heights). The final variances were = ±0.2 and <r p = ±20% 
which indicates that the two quantities are separable and observable. 

Covariance propagations were performed for each of the particle size 
ranges. The results of these propagations show numerically how the initially 
assumed variances of the aerosol parameters are reduced by the information 
gained in a measurement sequence. 
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Initial values of the aerosol parameter RMS uncertainties were chosen 
to reflect the range of values occurring in the literature (Newell and Gray, 
1972, and Malchow, 1971) # Number density (>°) was assumed to vary over 
one order of magnitude ±3cr, thus the a p used was 166%. This is consistent 
with _the large observed, variations in the number density of stratospheric 
aerosols related to volcanic activity. The size parameter (a) was assumed 
to cover the entire observed range 0<a<6, thus a Q = 1. The real part of 
the index of refraction ranges from that of water (1.33) to fused silica 
(1.65) thus yielding a ±3a range of 0. 32 and a a of 0.05. Finally the complex 
part of the index of refraction was assumed to be as large as 0.015 (+3 ct) 
which is felt to be conservative (Volz, 1973) with a . a , value of 0.005. 
Both the real and complex parts of the index of refraction are assumed to 
be wavelength independent. All initial RMS values are assumed constant 
with altitude. 

Figure 2.1-13 shows graphically the reduction in RMS uncertainty 
resulting from a covariance propagation for the (0.1- l.Oyu) particle size range. 
The solid line indicates the assumed initial value of a particular quantity, 
and the error bars show the final RMS uncertainty. Numerical values 
corresponding to the error bars are listed in Table 2.1-1. Considerable 
reduction of the number density uncertainty occurs for all altitudes although 
it is especially prominent at 10 and 20 km where the aerosol turbidity is 
at a maximum. A similar pattern of error reduction is displayed for the 
size parameter (a). Uncertainty in the complex part of the index of refraction 
is reduced, however, the improvement is not as striking as in the case of 
the other parameters. In general the covariance propagation for this size 
range shows that substantial reductions in the initial uncertainties of aerosol 
number density, imaginary index, and size parameters are effected by the 
inversion. The uncertainty in the real part of the index of refraction is 
however reduced only slightly. This result is consistent with the sensitivity 
curves in Figures 2.1-1 and 2.1-2 which show little sensitivity to the real 
index of refraction at any altitude. 


12 



2.1.3 (0.01 -0.1a/) Aerosol Inversion Simulations 

The aerosol number density was adjusted to preserve the standard 

extinction versus altitude curve for 5500 A with ce = 3, n = 1.5, and n' = 0. 

7 3 

For this range the density varied from 1.4x10 /cm at the ground to 5.5 x 
- 2 3 

10 , /cm at 100 km. Since the phase function for this size range varies 
less with angle than in the (0. 1-1.0//) range only a few zenith angles were 
run. The results of these runs show that the constituent sensitivities did 
not significantly change with zenith angle (see Figures 2.1-6 and 2.1-7). 

The most notable difference between this size range and the two others 
is the lack of sensitivity to (n 1 ). For these small particles only number 
density and the size distribution parameter (a) are observable because the 
partial derivatives of the scattering phase function with respect to (n) and 
(n') for this size range, are near zero. A change from 300% to 100% 
uncertainty in number density enhances the (a) sensitivity rather than 
increasing (n') sensitivity (see Figure 2.1-6 and 2.1-5 respectively). 

Analysis of the two sensitivity curves for ( fl ) and (a) indicates a 
potential problem in distinguishing between the two parameters. The fact 
that the sensitivities for both parameters have the same wavelength 
dependence implies that the measurement equations at different wavelengths 
are not totally linearly independent (see Figure 2.1-8). However a coarse 
covariance propagation was run (five wavelengths and four kilometer tangent 
height separation) which indicates that the two are separable. The 
correlations were indeed high (0.6 to 0.8) between the parameters but there 
is a factor of two decrease in number density variance and a factor of ten 
decrease in size distribution variance from the initial values. 

The two parameters ( P and a ) are highly observable from 10 to 100 

km or until instrument noise becomes dominant. Typically an instrument 
-3 ‘ 

noise of 10 times the maximum horizon signal restricts the altitude 

• -4 

sensitivity range to 85 km while a noise value of 10 extends the altitude 
sensitivity to approximately 100 km. Thus the interference that comes 


13 



about from increased uncertainties in neutral density at higher altitudes 
does not reduce the observability of the two aerosol parameters significantly. 

Figure 2.1-14 illustrates the RMS error reduction for the particles 
in this size range. The numerical values are listed in Table 2.1 -II." The 
results are generally similar to those obtained for the (0. 1-1. 0/0 size range 
except that there is a smaller reduction of the number density uncertainty. 
Somewhat surprising, in view of the near zero sensitivity for (n 1 ) shown in 
Figure 2.1-6, is the substantial reduction of the uncertainty in (n 1 ) after a 
covariance propagation. This indicates that the factor causing the error 
reductionisa buildup of correlations between (n 1 ) and the other parameters 
as opposed to a partial derivative effect which dominates the sensitivity 
results. When the sensitivity and covariance propagation results differ 
significantly as in this case it is important to consider the differences 
between the abbreviated covariance propagation used to produce these 
working numbers and a real covariance propagation used in an actual 
inversion. To insure the ultimate convergence of the state vector estimate 
it is often required to reduce the covariance update by means of a numerical 
multiplier. This results in a smaller reduction of the initial covariance 
after a given number of reiterative updates than would be indicated by the 
unweighted propagation used here. 

2.1.4 (1.0-10//) Aerosol Inversion Simulations 

The large particles in this range are the furthest from Rayleigh in 
their scattering properties, i.e., in their cross-section wavelength 
dependence and in the shape of their phase function. Because of their large 
cross section the number density is low, running from 23.4 /cm at the 
ground to 7.7x10 / cm^ at 100 km. 

The great difference between these particles and other atmospheric 
constituents makes them the most readily inverted of all size ranges. It is 
possible to invert all four aerosol parameters from a given horizon scan 
when the phase function is not minimized by an unfavorable scattering angle 
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i.e. sun angle. Although the four sensitivities are reasonably close, results 
from the previous size ranges indicate that all four physical characteristics 
are separable. Variable scattering singles are readily obtainable for nearly 
all conceivable orbits from equatorial to polar orbits with the exception of 
the 600 and 1800 local hour angle near polar-orbits. 

There is, however, a problem with this size range that does not appear 
in the other ranges. That is, when the scattering angle approaches 110°, 
the phase function is so small that the aerosol energy contribution and 
respective sensitivities drop drastically and only number desity is highly 
invertible from the data (Figure 2.1-12). The variations in sensitivity with 
changing solar angle can be clearly seen in Figures 2.1-9 through 2.1-12. 

Covariance propagation results for the large particles are shown in 
Figure 2.1-15 and the corresponding numerical values are listed in Table 
2.1 -II. In the case of the large particles all the parameter uncertainties 
are substantially reduced including the real part of the index of refraction 
which was essentially unaffected for the other particle size ranges. 

2.1.5 Modeling Errors 

With the standard filter approach to the inversion problem it is 
theoretically possible to drive the variances on all the parameters as close 
to zero as desired by taking enough measurements. However, since there 
are inherent errors in the radiative transfer modeling a near perfect variance 
is unrealistic. The purpose of this section is to determine the effects of 
aerosol modeling errors on the final accuracy obtainable in an inversion. 
This is done by computing the vectors (k, Ax) and the scalar (A I) in the 
equation Ax = k (AI) where Ax is the state error, (k) is the filter gain, and 
AI is the intensity error. 

The results are presented in sets of tables. The first table of each 
set (Table 2.1-IV) illustrates the gain (k) which is used to compute errors 
in the state given errors in intensity. This is done with a geometry having 
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a zenith angle = 30° and an azimuth angle = 0° for all size ranges, with 
standard initial errors (see section 2.1.2), and using a 300% initial uncertainty 
in aerosol number density. The additional tables of each set are the computed 
state error Ax for errors in the albedo and aerosol models and for an 
instrument bias. The state errors are based upon the gain (k) of the first 
table. 

It is important to realize that the initial uncertainties are important 
in the propagation of error. For example, a parameter which has both a 
large partial derivative and a large initial uncertainty will receive a large 
update and become more in error by inaccurate models than a parameter 
with a small initial variance and partial derivative. The numbers given 
represent the maximum error that will occur in an inversion using 3000, 
4000, 5500, and 7000 A as wavelength channels. It should be noted that the 
(k) values for 3000 A at 20 and 40 km, and for 4000 A at 20 km are set 
equal to zero. This is because these wavelengths are saturated at the 
designated altitudes and therefore yield no information about the densities. 

Tables 2.1 -V, VI, IX, X, XV, and XVI list the errors in the state 
vector elements produced by errors in the aerosol extinction and phase 
function models. The model errors used here are representative of the 
errors associated with the aerosol models presented in Section 3.3 of this 
report. Tables2.1 -V, VI list the state errors for aerosols in the (0.01 -0.1 a/) 
size range. The modeling errors for this size range are small («2% in 
extinction and « 5% in phase function) and consequently the effects on the 
state are also small. The largest resulting error is due to the phase function 
error, and is a 12.3% error in the particle number density at 60 km as 
shown in the fourth column of Table 2„1 -VI. 

The (0.1-I.Oa/) particle size range has larger associated errors and 
consequently larger effects on the state. Tables 2.1 -IX, X list the state 
errors associated with this size range. As with the smaller size range, 
the phase function error (10%) produces larger state errors than the 
extinction error (5%). The largest error induced is in the aerosol number 
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density at 60 km, and is 26%. This error is still small compared to the 
assumed initial number density uncertainty of 300%, and is less than the 
60% uncertainty obtained after a coarse covariance propagation. 

For the large particle range the state errors are presented in Tables 
2.1 -XV, XVI. Again the number density estimates are primarily affected, 
and the errors are comparable to the final state estimates obtained via a 
covariance propagation. 

Table 2.1 -XII shows the effects of the assumption of the wrong size 
range for aerosols, i.e., the table shows what errors result if the size 
range is assumed to be ( 0.1-1. 0/u) whereas the actual size range is (,01-0.1/u). 
The resulting errors are large as expected since the extinctions differ by 
two orders of magnitude for. the two size ranges. Aerosol number density 
is shown to be in error- by one order of magnitude or « 1000% at 40 km. 
Other state elements are also strongly affected. The neutral density estimate 
is in error by 52% at 80 km, and the ozone density by 242% at 60 km. 
Errors in the other aerosol parameters are significant though not large 
(these are given in absolute units). At 40 km, the error in (a) approaches 
one which represents a significant alteration of the size distribution. The 
results of this size range shift emphasize the importance of either knowing 
rather accurately the particle size limits, or expanding the state vector to 
include these limits. 

Table 2.1 -XIII shows the effects of the introduction of a 10% instrument 
bias error in each wavelength channel. The error induced in the aerosol 
number density is moderately large because the signal contribution by the 
aerosol to the total signal is small. Since the aerosol number density 
uncertainty is large, the filter attempts to correct the signal error by 
adjusting mainly the aerosol number density. 

Finally, Tables 2.1 -VII, XI, and XVII show the effects on the state 
of the introduction of an unestimated effective surface albedo deviation from 
the expected value. The effect of an albedo uncertainty is similar to the 
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effect of an instrument bias, i.e., there is a broad-band shift in the measured 
intensity away from the expected value. Since aerosol density is assumed 
to have a relatively large initial uncertainty, it is this quantity which is 
mainly adjusted to account for the intensity shift. The amount of adjustment 
done to the aerosol number density is inversely proportional to the 
contribution that the aerosol makes to the total intensity. Error values in 
Table 2.1 -VII, XI, and XVII are based upon a scattering angle that coincides 
with the minimum in the angular scattering functions for the medium and 
large-sized particles. 

Thus these aerosols are making a minimal contribution to the total 
intensity at the receiver, and therefore the errors related to an uncertain 
albedo are maximized. Tables 2.1 -VII, XI, and XVII illustrate this effect 
dramatically. Each aerosol size range is defined to have the same extinction. 
As the particle size increases, the scattering phase function decreases (at 
the chosen sun angle) and the error increases rapidly for the same extinction 
in each size range. There are two implications of these results. One is 
that the albedo should be added to the state vector and estimated to minimize 
its uncertainty. The other implication is that the solar scattering angle 
for primary radiation should be chosen to minimize the effect of albedo 
uncertainties. At more desirable scattering angles the large particle errors 
(Table 2.1 -XVII) are substantially reduced, and are comparable to the 
small particle errors (Table 2.1 -VII). For example, the 1160% error in 
number density for large particles at a scattering angle of 80° (Table 2.1 
-XVII, 40 km) is reduced to 60% when the scattering angle is changed to 
30°. 
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TABLE 2.1 -IX 
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Fig. 2.1-1 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Range (0.1-1.0y), Zenith Angl 
Azimuth Angle = 0°, and la values (p=100%, a-1 


Wavelenath in Anostrnmc 



Fig. 2.1-2 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Ranqe (O.l-l.Ou), Zenith Angle = 60°, 
Azimuth Angle = 0° , and la values (p=300%, a=l , n=0.05 
n' =0. 005 ), Wavelength in Angstroms. 
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Fig. 2.1-3 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Range (O.l-l.Oy) , Zenith Angle 
Azimuth Angle = 0°, and lo values (p=600%, a=l, 
n'^.005), Wavelength in Angstroms. 
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Fig. 2.1-4 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Range (O.l-l.Oy), Zenith Anale 
Azimuth Angle = 0°, and lo values (p=300%, a=l, 
n'=.005), Wavelength in Angstroms. 
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Fig. 2.1-6 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Range (0.01-0. Ip) , Zenith Angle = 6 
Azimuth Angle = 0°, and lc values (p=300%, a=l , n=.0 
n'=.005), Wavelength in Angstroms. 
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Fig. 2.1-9 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Range (1.0-lOy), Zenith Angle = 60°, 
Azimuth Angle = 0°, and la values (p=300%, <S=1, n=.05) 
n'=.005), Wavelength in Angstroms. 
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Fig. 2.1-11 Scattered Sunlight Constituent Sensitivities 

for Aerosol Size Range (1.0-10u), Zenith Angle 
Azimuth Angle = 180°, and la values (p=300%, a 
n=.05, n’=.005), Wavelength in Angstroms. 
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Fig. 2.1-13 Covariance Propagation for Aerosol Size Range 

(0.1-1. Op), Zenith Angle = 60°, Azimuth Angle = 
and lo values (p=166%, a=l, n= . 05 ,n ' =. 005) , Wave 
lengths of 3000, 4000, 5000, 6000, and 7000 8, 
and Tangent Heights Spaced at 4 km. 
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Fig. 2.1-14 Covariance Propagation for Aerosol Size Range 

(0,01-0.1y) Zenith Angle = 60°, Azimuth Angle = 0 
and la values (p=166%, a=l, n=.05, n'=.005), Wave 
lengths of 3000, 4000, 5000, 6000, and 7000 X, 
and Tangent Heights Spaced at 4 km. 
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Fig. 2.1-15 Covariance Propagation for Aerosol Size Range 

(T. 0-10y) , Zenith Angle = 60°, Azimuth Angle = 0° 
and lo values (p=166%, a=l , n=.05, n'=.005). Wave 
lengths of 3000, 4000, 5000, 6000, and 7000 &, 
and Tangent Heights Spaced at 4 km. 



2.2 STELLAR OCCULTATION SIMULATION 

2.2.1 Introduction 

A One of the contract- subtasks was concerned with the application of 
the Kalman- Bucy filter to the inversion of stellar occultation data. The 
data, however, was unobtainable during the contract period. The inversion 
scheme was nevertheless developed and a sensitivity analysis was 
performed. Several computer runs were made using simulated data which 
illustrates the invertibility of high altitude aerosol layers with a star 
occultation experiment. 

2.2.2 Sensitivity 

The sensitivity function for stellar occultation is identical to that 
used for scattered light. The partial derivatives are, however, applied to a 
different radiative transfer model. Figure 2.2-1 shows the envelope of 
peak sensitivities for a complete stellar occultation scan at several 
wavelengths. It can be seen from this figure that both ozone and neutral 
density will be invertible down to an altitude of 20 km. The curves indicate 
ozone cannot be inverted below this level. This is because the intensity 
attenuation of ozone sensitive wavelengths below 20 km is too great to permit 
measurement of the signal with the assumed instrument signal to noise 
ratio of 100:1. 

Figure 2.2-2 shows the effect of a high altitude aerosol layer (50 
km, one order of magnitude anomalous increase) on the stellar occultation 
sensitivity curves. It can be seen that standard levels of aerosols do generate 
a sufficient signal increase to be. detected, however, an anomalous increase 
of thisnature will be measured and inverted for the assumed 100:1 signal 
to noise ratio. 
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2.2.3 Selected Simulation Results 


For the purpose of simulation, fictitious, but realistic, density 
distributions other than standard were used to generate simulated data e.g. 
anomalous layers were included based on measurements of Rossler, 1968 
and Elliot, 1971. The data was then processed by the inversion routine. 
In Figure 2.2-3 the crosses represent the estimate of the density with two 
error bars given for each point. The horizontal bar represents the un- 
certainty in the altitude at which the average density occurs while the 
vertical bar represents the density uncertainty. The density error 
bars are derived from the filter while the altitude error bars are derived 
a priori in simulations. The solid line represents the fictitious density or 
"right" answer. If the densities fall within the RMS error bars 67% of the 
time the simulation is considered a success. 

Figure 2.2-3 is a simulation in which an order of magnitude anomalous 
aerosol layer at 50 km was used to generate fictitious data. The inversion 
results demonstrate the ability of the inversion technique to determine such 
layers if they exist with a 100:1 signal to noise ratio and wavelengths of 
3000, 4000, 5000, and 7000 A. 

2.2.4 Conclusion 

The Kalman-Bucy filter is easily adapted to the inversion of stellar 
occultation data and produces excellent results. Although a simplified 
radiative transfer model was used in these simulations there is no reason 
why a more sophisticated model incorporating refraction and dispersion 
effects could not be used. This would extend the range of invertibility well 
below 30 km for both neutral density and aerosols in addition to providing 
more accurate results above 30 km. 
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SENSITIVITY DENSITY PEAK 


7000 



Fig 2.2-1 Stellar Occultation Constituent Peak Sensitivity 
Profiles with 100:1 signal to noise and standard 
densities . 
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7000 



Fig. 2.2-2 Stellar Occultation Constituent Peak Sensitivity 
Profiles with 100:1 signal to noise and an (order 
of magnitude increase) anomalous aerosol layer at 
50 km. 
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Inversion of Simulated Stellar Occultation Data in 
the presence of an anomalous aerosol layer at 50 km 





2.3 HORIZON PROFILE PARAMETER SIMULATION 
2.3.1 Introduction 

In this section we discuss the results of a horizon profile parameter 
study aimed at determining how horizon profiles are altered by changes in 
the various measurement parameters. These include; season, latitude, 
satellite altitude, solar zenith and azimuth angles, and either the earth's 
albedo or the altitude and albedo of a cloud layer. The newly developed 
mulitple scattering code (C.K.W. ) was used to produce 324 radiance profiles 
representing various combinations of the measurement parameters. The 
horizon profiles are plotted in multiples of from three to five profiles per 
graph and are compiled in Appendix 5.1. Figure 2.3 _ 1 is a sample profile 
set showing zenith angle as a variable. The measurement parameter being 
varied in each graph is identified at the top of the graph and values of a 
given variable corresponding to a given profile are listed to the right of 
each graph. The fixed parameters for each profile set are listed in the 
legend. The right hand graph in each figure displays the fractional change 
of each radiance profile relative to one member of . the set. 

Figure 2.3- 2 shows the limb scan geometry, defining the tangent height 
(H), the satellite altitude (S), and. the line of sight (L). The direction of 
the sun's rays are specified in relation to a Cartesian coordinate frame 

A A 

attached to the satellite, with (Z) along the local vertical and with (Y) in 

A A 

the plane defined by (Z) and (L). The solar zenith angle (0) is measured 

A A A 

relative to (Z) and the solar azimuth angle ( <f> ) is measured in the (X,Y) 

A A A 

plane relative to (Y). The angle which (L) makes with (Y), i.e., the scan 
angle, is defined as (6). 

All of the horizon profiles given in this report pertain to a satellite 
altitude of 500 km. The scan angle (5) therefore varies from about 19.5° 
to 2 2° as the tangent height varies from 1 00 km to zero. The single scattering 
angle is determined by the angles (0, <f>, and 6 ) through the expression cos<^ 
= (cos0sin0cosd - cos0sin<5). Thus, for the coplanar cases (<£=0° and <^ = 180°) 
0=90 °+S~e and 0=9O°+<5+0, and for the <f>= 90° case cos^ is simply given by 
(-cos0sin6). 
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2.3.2 Discussion of Results 


For reasons of organizational convenience the 84 figures for this 
parameter study are located in Section 5.1. In each case the atmosphere 
was considered to have two-kilometer layers (DZ = 2) that extended to an 
altitude of one hundred and twenty two kilometers (ATM = 122). The seasonal 
and latitudinal distributions for the mean neutral and ozone densities used 
in determining these profiles are shown in Figures 5. 1 - 73 to 5.1-84. Figures 
5.1-1 to 5.1-12 show the variations in horizon profiles resulting from 
seasonal variations in the neutral and ozone densities for each of twelve 
latitide- wavelength conditions (latitude = 0, 40, and 70°N and \= 3000, 
4000, 5500 and 7000 Angstroms). At 0°N latitude the same neutral density 
distribution was used for all seasons. The changes in the profiles with 
season at 0°N latitude reflect, therefore, only the changes in the ozone 
density distribution. This profile set represents the only case where the 
differential horizon profiles can be interpreted simply in terms of the 
changes occurring in iust one of the atmospheric constituents . These changes, 
which are referenced to the winter profile, are greatest at 3000 A and 
smallest at 4000 A in agreement with the known spectral characteristics 
of the ozone absorption cross section (Wu, 1970). The remaining horizon 
profiles in this parameter study are arbitrarily based upon the latitude 
and seasonal choices of 40°N latitude, winter. 

Figures 5.1-13 to 5.1-24 show the variations in horizon profiles with 

wavelength for twelve cloud- zenith angle conditions (a zenith angle of 30° 

and 80° and a cloud at 4, 7, and 10 km with an albedo of 0.4 and 0.8). The 

solar irradiances at 3000, 4000, 5500, and 7000 A are 5.14, 14.92, 17.25, 

2 

and 13.25 /vwatts/cm -A respectively (Thekaekara and Drummond, 1970). 
Figures 5.1-25 to 5.1-28 show the variations in horizon profiles with 
wavelength for four zenith angles (0, 30, 50, and 80 degrees) and a ground 
albedo of 0.3. These figures show that clouds are only observable at the 
longer wavelengths like 5500 and 7000 A which penetrate deepest into the 
atmosphere. In addition, these figures show that the enhancement of radiance 
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due to an increase in cloud albedo is strongly dependent upon the zenith 
and azimuth angles, which determine the scattering angle (0) and hence 
the fraction (sin0) of the solar irradiance that is actually incident upon the 
cloud surface. 

Figures 5.1-29 to 5.1-64 show the variations in horizon profiles with 
zenith angle for 36 combinations of azimuth angle, wavelength, and ground 
albedo (azimuth angle = 0, 90, and 180 degrees; A= 3000, 4000, 5500, and 
7000 Angstroms; ground albedo = 0.2, 0.5 and 0.8). The differential horizon 
profiles are all referenced to the zero zenith angle profile. These figures 
show that the shape of a profile is not affected greatly by changes in the 
albedo and in the solar azimuth and zenith angles. For the brightest profile 
(80° zenith and 0° azimuth) there is little albedo effect because sin (0) is 
relatively small and therefore most of the signal is produced by strong 
forward scattering high in the atmosphere. Albedo effects are most 
prominent for azimuth and zenith angle combinations that make the scattering 
angle (0) ~ n/2. 

Figures 5.1-65 to 5.1-72 show the variations in horizon profiles with 
ground albedo for eight combinations of zenith angle and wavelength (zenith 
angle = 30 and 80 degrees; X= 3000, 4000, 5500, and 7000 Angstroms). 
The differential horizon profiles are referenced to the albedo = 0 profile. 
These figures show that the enhancement of the line of sight radiance due 
to the earth's albedo is strongly dependent upon wavelength and scattering 
angle. 


51 






Fig. 2.3-2 Limb Scan Geometry 
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3.0 


MODELING 


3.1 FILTER DEVELOPMENT 

3.1.1 Introduction 

— y — ... 

There are a number of approaches to the problem of inverting aerosol 
parameters, neutral atmospheric density, and ozone density from horizon 
profile data. Ideally the method chosen should provide a best estimate (in 
some statistical sense) of the desired parameters incorporating not only 
the data, but the instrument characteristics arid previous knowledge of the 
parameters as well. What follows is a discussion of several techniques 
which have been considered as well as the reasoning behind the choice of 
the Kalman- Bucy filter. 

3.1.2 Comparison of Techniques 

The first and most obvious approach to the inversion problem is to 
attempt an analytic solution of the equations of radiative transfer, that is 
solving for density as a function of intensity. This technique, has been 
used (Anderson, 1969) for atmospheric probing; however, even for simple 
geometries with single scattering and only one unknown constituent several 
simplifying assumptions must be made in order to reach an analytic solution. 
Also, once the solution is reached it is difficult to perform a meaningful 
error analysis on the final answer. Thus given the stated goals of the 
method to be chosen and the complexity of the radiative transfer equations 
for horizon scan geometry, it is clear that this approach to the inversion 
problem is limited. 

Since it is not practically feasible to solve the equations analytically 
the next possibility to consider is an indirect solution to the equations using 
partial derivatives and some iterative scheme to arrive at an approximate 
solution to a set of simultaneous equations. This would probably be possible 
and would give reasonable answers, however, there would be no statistical 
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error analysis. Also if there were more data points than needed to solve 
the simultaneous equations only the first data taken would be used. Thus 
this technique involves incomplete data utilization and does not provide 
for the desired error analysis. 

Next, a regression technique might be considered. A regression would 
provide not only a best estimate in a least squares sense but error statistics 
as well. However there are some problems with the application of a 
regression and these are: 1) Previous knowledge of the parameters is not 
incorporated in the estimate or the statistics, 2) Instrument noise is not 
incorporated, and 3) The large amount of data must be processed 
simultaneously rather than sequentially. However this would be a good 
technique if there was not a better approach to the problem. 

The use of Bayesian statistics makes possible the incorporation of 
previous knowledge and instrument noise in the inversion scheme. It also 
allows data to be processed one point at a time thus incorporating each 
new datum into the estimate of the statistics of the estimate. There are a 
number of Bayesian approaches to the problem such as maximum likelihood 
estimation, Bayesian regression, weighted least squares, and the Kalman- 
Bucy filter. 

For the case in which the equations relating the measurement and 
the state are linear, and noises are Gaussian, all these techniques are 
equivalent. However in real physical situations these assumptions are often 
invalid and the techniques differ. The Kalman- Bucy filter was chosen 
because: 1) It has all the advantages of Bayesian statistics, 2) The filter 
equations are simple and easy to understand and code, and 3) Because of 
the simplicity of the equations the filter readily lends itself to iterative 
techniques which ensure convergence in spite of the nonlinearity of the 
equations. 
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3.1.3 Filter Review 


The Kalman- Bucy filter is a recursive technique for estimating the 
state (in this case aerosol parameters, neutral atmospheric density, and 
ozone) of a system (the atmosphere) from measurements of that system 
(horizon profile data). At each recursion the previous estimate of the state 
is updated. The degree and distribution among parameters of the update 
is determined by the measurement the previous estimate of the state, the 
instrument noise, and the statistics (covariance matrix) associated with the 
estimate of the state. After the state is updated the covariance matrix is 
also updated so that it reflects the current knowledge of the state. 

The updating procedure is done in three steps. First is the calculation 
of the filter gain which determines the degree and distribution of the update 
among the elements of the state vector. 

K(m+1 ) = P(m) B T (m+l) [B(m+1) P(m) B T (m+l) + R]' 1 (3.1-1) 

where 

(m) is an index of measurement 

(K) is the filter gain 

(P) is the covariance matrix 

(B) is the measurement vector whose elements are 
b^ = 3h(x, m+l)/3x^ 

h(x, m+1) is the predicted value of the intensity based on the 
state estimate and measurement parameters 

and 


(R) is the instrument noise. 
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The second step is the updating of the state estimate 


x(m+l) = x(m) + K(m+1) [i (m+1) - h(x, m+1)] (3.1-2) 

where (x) is the state estimate and (I) is the measurement. 

Thirdly the covariance matrix is updated to incorporate the knowledge 
gained of the state. 

P(m+1) = P(m) - K(m+1) B(m+1) P(m) (3.1-3) 

As can be seen from the equations there are three steps which must 
betaken before applying a filter to the problem of limb scan data inversion. 
First a computer code giving a direct solution to the equations of radiative 
transfer is needed for the calculation of h(x, m) and 3h(x, m)/3x (partials 
are usually calculated by the approximation Ah(x, m)/Ax). Second the state 
(x) must be defined, and third some initial estimate of the covariance matrix 
P(0) must be made to start the recursive process. 

Two computer codes are available for the first requirement and these 
are discussed in detail in Section 3.2. For all of the inversion work to 
date the hybrid single scattering code has been used since the multiple 
scattering code had not been completed. However, the characteristic shape 
of the hybrid profile agrees well with the multiple scattering code. Also, 
the hybrid codes use less computer time and it is simpler to extract the 
partial derivatives from the single scattering code. It is important to note 
that one of the advantages of the filter inversion technique is the capability 
of incorporating any solution into the equations of radiative transfer no 
matter how complex. Thus in an actual experiment the more accurate 
multiple scattering code would be used. 

The definition of the state has been previously discussed (Newell and 
Gray, 19 7^) for the purposes of inverting three constituent densities. 
Basically that scheme had density points defined at arbitrary altitude 
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increments. The elements of the state were the logs of the densities (or 
extinction in the case of aerosols) at each altitude. By using this definition 
of the state at intervals of from 3-10 km along with a finer integration 
step of 1 km in the radiative transfer codes it was found that the filter 
would converge (Newell and Gray, 1972) for realistic densities that would 
be encountered in the atmosphere. 

3.1.4 Sensitivity Review 

For this current work the state definition was expanded to include 
the log of aerosol number density ( P ), size distribution parameter (a), and 
both real and complex indices of refraction (n and n 1 ) at each chosen altitude. 
This expanded state was first analyzed through sensitivity studies and then 
through covariance propagation to determine the feasibility of estimating 
each parameter. The results of this study are discussed in Section 2.1. . 


The sensitivity vector is derived from the filter equations and indicates 

the relative information derived about each constituent parameter at each 

altitude for every measurement condition (i.e. wavelength and viewing angle). 

Since after each measurement the covariance matrix records the information 

gained from the measurement it is natural to look at the covariance update 

equation for the definition of sensitivity. The term of interest here is 

K(m+1) B(m+1) P(m), which gives the decrease in the covariance after the 

(m+1) measurement. Normalizing by P(m) yields the sensitivity matrix 

K(m+1) B(m+1) whose diagonal elements comprise the sensitivity vector. 

For practical data presentation and speed in calculation, the sensitivity 

vector can be defined by assuming that the covariance matrix maintains 

its initial variances and has zero cross correlation. With these assumptions 

th 

the sensitivity of the j n element of the state vector is 


Sj (m) 



/ 3h(x,m ) \ 2 

V 8x i / 



+ R 


(3.1-4) 



It is important at this time to realize both the usefulness and drawbacks 
of sensitivity analysis i The main purpose of the analysis is to show just 
which constituent parameters and altitude regions are important in a 
measurement. Thus sensitivity analysis will reveal for example that at 
4000 A and a tangent height of 20 km with a certain set of initial uncertainties 
the greatest information will be derived about aerosol number density and 
aerosol absorption in the region from 20 to 23 km. The analysis will also 
reveal that at 3000 A and 20 km the only information derived is from ozone 
between 40 and 50 km. (This is due to the strong ozone absorption which 
prohibits this wavelength from penetrating below 30 km.) 

The one disadvantage of the sensitivity analysis stems from the 
assumption that each update has zero' cross correlation. It often happens 
that a set of measurements do not yield linearly independent sensitivity 
vectors. For example consider a simple experiment designed to determine 
aerosol number density and neutral density at 50 km by two measurements 
at 3800 and 4000 A. The sensitivity vectors would be high for both 
constituents at both wavelengths seemingly indicating that both are readily 
observable. However, the actual final covariance matrix of such an 
experiment would show that little has been learned of either constituent 
since there is no way to distinguish between constituents from the 
measurements. Therefore, the covariance matrix would indicate that the 
sum of both constituents is well known. 

, This is not a serious problem however since visual inspection of 
sensitivity vectors is often sufficient to determine whether or not there 
will be a problem with ambiguous measurements. If there seems to be 
some doubt then a complete propagation of the covariance matrix can be 
made to determine exactly what information is recoverable from the data. 

3.1.5 Stellar Occultation 

The stellar occultation inversion technique is essentially the same 
as the inversion of scattered sunlight (Newell and Gray, 1972). There are, 
however, some aspects of the stellar occultation radiative transfer model 
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that make improvements possible in the construct of the state elements. 
These improvements yield finer data resolution and more accurate error 
statistics. The next few paragraphs describe the essential differences 
between the two inversions. 

. The radiative transfer model used for the stellar occultation work 
is 

1 = I o e" T (3.1-5) 

where (I) is the measured intensity, (I Q ) is the star intensity, and (r) is 
the optical depth along the line of sight. The model used assumes a linear 
line of sight and neglects dispersion and multiple scattering effects., These 
assumptions will not cause errors in excess of 10% above 30 km (Hayes 
and Roble, 1968). This altitude is approximately the maximum depth to 
which ozone information can be retrieved. ■ , ■ ■ 

Since the Kalman-Bucy filter assumes linearity in the measurement 
equation it is desirable to obtain a radiative transfer model which is as' 
nearly linear as possible. In the case of stellar occultation it is possible 
to, create an exact linear equation by measuring the log of the intensity 
rather than the intensity thus giving the new measurement equation 

ln(I) = ln(I Q ) - r (3.1-6) 

In the scattered light case discussed previously it is necessary, to 
estimate densities at a resolution approximately two to three times as coarse 
as the tangent height data. This was necessary along with a gain damping 
term to insure convergence of the filter. For the linearized stellar, 
occultation case it is -possible to estimate density at each tangent height 
and also to dispense with the gain damping factor. The final density estimates 
are therefore average densities of a layer between two tangent heights. 
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3.2 RADIATIVE TRANSFER MODELING 

3.2.1 Introduction 

This section discusses the extension of computer algorithms for 
simulating the horizon profile caused by radiative transfer within the earth's 
atmosphere (previously reported in Newell and Gray, 1972). Here, we begin 
by reviewing the computational problems to which these codes are addressed 
and introduce briefly computer codes applicable to related problems. Their 
applicability to present work objectives is discussed, with emphasis on 
the unsolved problems that required the development of new codes. These 
new codes, a simple, modified single scattering code and a more 
sophisticated, full multiple scattering code, are compared in terms of 
accuracy obtainable versus computational burden required. The various 
applications to which they have been put are presented. 

3.2.2 The Computational Problem 

The problem of horizon profile simulation is a difficult one because 
in the earth's atmosphere there are important polarization effects in the 
presence of multiple scattering involving not only Rayleigh but also aerosol 
scatterers, as well as absorption. Furthermore, there is the 
nonhomogeneous structure of the atmosphere, and its slight curvature in 
conformity with the earth. Many of these factors contribute substantial 
difficulty. Polarization requires a matrix rather than a scalar treatment, 
thus increasing the complexity of operations and the amount of storage 
space required. Multiple scattering is a well-known primary source of 
difficulty in all radiative transfer calculations. Rayleigh* scattering and 
absorption are straightforward, but aerosol scattering is characterized by 
an exceedingly complicated ill-behaved phase function, or angular pattern, 
that is difficult to incorporate in a radiative transfer model. The 
nonhomogeneous structure of the atmosphere imposes an altitude dependence 
of the scattering into any given solid angle. Finally, the curvature of the 
earth entails formidable geometry problems. 
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3.2.3 Codes Applicable to Related Problems 


Basic computational techniques for radiative transfer modeling that 
existed at the outset of our work are reviewed by Hunt (1971). We shall 
refer to the two basic approaches he discusses as the matrix operator 
method and the Monte Carlo method. The matrix operator method represents 
the atmosphere (or any sub- layer of it) as a matrix operator that transforms 
an input column vector of stream irradiances to an output column vector 
of stream irradiances.- This matrix formulation is not to be confused with 
the matrix formulation required by the inclusion of polarization in the 
problem; the matrix operator method has matrices even if polarization is 
ignored. The Monte Carlo.method represents the atmosphere with a random 
number generator that chooses scattering locations, angles etc. for individual 
photons, very large numbers of which are followed to accumulate a horizon 
profile. - 


The matrix operator method has been used mainly in homogeneous 
plane parallel atmospheres, where application of the doubling technique 
makes it impressively efficient. It has not, however, yet been developed 
to handle inhomogeneous curved atmospheres. Thus, while it was available 
in principle at the outset of our work, in practice, its use would have required 
extensive development. By contrast, the Monte Carlo method has actually 
been applied to essentially our problem in a program called FLASH written 
by Collins and Wells, 1970. The FLASH program will indeed accommodate 
every one of the sources of difficulty mentioned in the preceeding section. 
Furthermore, the program existed at the outset of our activities. It would 
have been an attractive solution to our problems except for one factor: it 
requires a large amount of computing time (approximately 2000 seconds 
of CDC 6600 cpu time per horizon scan to achieve an accuracy of 
approximately four percent). 
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3.2.4 The Radiative Transfer Codes 


Within the Draper Laboratory, there are now two radiative transfer 
codes that simulate horizon profiles. One is a refinement and completion 
of the code REV. HYBRID which is described briefly in the earlier report 
by Newell and Gray, 1972 and in greater detail in an internal report by 
Var, 1971. It is based essentially on single scattering plus an effective 
planetary or cloud albedo as a secondary source. The second code (C.K.W.) 
is a more recently developed full multiple scattering model with polarization. 
It is based on mathematical techniques described in Newell and Gray, 1972 
and in more detail in a recent journal article by Whitney, 19 72. Some of 
the computational procedures involved in implementing the model are 
described in Appendix 5.3 of this document. 

Considerable effort has been expended to provide a basis for judging 
which of the two radiative transfer codes should be used for any particular 
application. The remainder of this section reports on these activities. 
There are three main areas of comparison; accuracy achievable, run time 
requirements, and suitability for specific applications. 

In attempting to establish the absolute accuracy of the C.K.W. code 
we compared its results to FLASH results through the courtesy of the Air 
Force Cambridge Research Laboratory. Differences between FLASH and 
C.K.W, were found to be small and related primarily to fluctuations in the 
Monte Carlo results or differences between the details of layer definitions 
in the two codes. 

The run times for both the C.K.W and the REV. HYBRID codes appear 
to be significantly less than those required for Monte Carlo simulations, 
in spite of the fact that specific comparisons are complicated by a number 
of operational and statistical factors. Our comparisons have been based 
rather arbitrarily on the times required to produce horizon profiles which 
are in apparent agreement with each other. Typically, the C.K.W code 
produces a horizon profile in approximately twenty seconds, and 
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RE V. HYBRID in. approximately seven seconds, both on an IBM 3 60/7 5y while 
the Monte Carlo, code of FLASH uses approximately two thousand /'seconds 
of CDC 6600 time to produce comparable results • The significant comparison 
here, is between the C.K.W. code and the Monte Carlo code, since only 
they are really addressed to the full multiple scattering problem. A tentative 
estimate of : the time saving with C.K.W; code is a factor of Several hundred. 


. In general, one would anticipate that multiple scattering would increase 
observed radiances ’to those that have 1 been scattered directly into the 
■.receiver field of : view. The magnitude of that contribution as determined 
by the C.K.W. r code is- illustrated in Figures 3; 2- 1 and 3.2-2. The 
enhancement of the radiance occurs without substantially altering the profile 
shape. Such an effect can be simulated also with the REV. HYBRID code 
by introducing a fictitious planetary albedo. The ficticious albedo required 
to match multiple scattering results depends upon sun angle and wavelength. 
Figures 3.2-3 and 3.2-4 show the effects of an albedo range of zero to 
unity upon REV. HYBRID horizon profile at 4000 and 5500 A. At 5500 A 
the vertical optical depth of the atmosphere is smaller than at 4000 A, 
thus more sunlight reaches the surface and is reflected back into the 
atmosphere. The reflected radiation is scattered into the receiver field 
of view adding substantially more radiation than when the albedo is zero 
causing a significant variation in the enhancement for the two wavelengths. 

There is no way that the enhancement required to make REV. HYBRID 
conform to the C.K.W. code can be derived from basic principles. Thus 
we feel it is wise to use the full multiple scattering code (C.K.W.) for 
applications requiring absolute radiances, rather than simply horizon profile 
shape. Where only shape is required, REV. HYBRID is sufficiently accurate. 

We come now to the question of suitability for different applications. 
There are two important applications to consider here: the generation of 
horizon profiles for the inversion simulations used to produce the aerosol 
invertibility results of Section 2.1, and the generation of the 324 sample 
horizon profiles in Appendix 5. 1 illustrating the effects of varying wavelength, 
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latitude, season, cloud height, zenith and azimuth angles. For the inversion 
simulations, the single scattering REV. HYBRID was judged suitable because 
the simulation does not require absolute radiances but rather fractional 
changes in radiances induced by small perturbations in each of the quantities 
defining elements of the state vector. Furthermore, the calculation of these 
differential changes would have required the calculation of many complete 
horizon profiles, and hence a considerable amount of computer run time 
with C.K.W. This is not the case with REV. HYBRID; differential changes 
can be obtained without recalculating complete profiles. By contrast, for 
the parametric horizon profile study, the computer code (C.K.W.) was used 
because the effects of the parameters to be illustrated could not have been 
accurately obtained from REV. HYBRID. 
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3.3 AEROSOL MODELING 


3.3.1 Introduction 

The application of recursive filter algorithms to limb profile inversion 
of aerosols requires that at each recursion the aerosol optical properties 
be recomputed given an updated physical characteristic estimate. Aerosol 
optical and physical properties are related theoretically through the Mie 
equations which express the angular scattering contributions and cross 
section as a function of the physical parameters such as size and index of 
refraction. At each recursion the filter requires recomputed values of 
these quantities and also the partial derivatives of radiant intensity with 
respect to those physical characteristics that are included in the state vector. 
For example the partial derivative of intensity with respect to index of 
refraction is of the general form 3l/3n = (3I/3t) (3r/3n) + (3I/3P ($)) 

(3P (0)/3n). Evaluation of 3r/3n and 3P (0)/3n involves the use of detailed 
aerosol models. Since computational times on the order of a minute are 
involved in determining these quantities using the Mie series, the total 
computer time required for inverting a single scan would be on the order 
of hours per scan just for aerosol calculations. Such an amount of time 
would be out of proportion to other calculational requirements posed by 
the filter algorithm, and inconsistent with the goal of near- real-time 
inversion capability. 

The implications of this calculational complexity when matched with 
the desire for computational economy are that multicalculations with the 
Mie series are to be avoided, and some alternate procedure must be provided. 

3.3.2 Choice of Models 

Several alternate procedures for the Mie computations are suggested. 
For the simple one parameter size distribution with which we are concerned, 
there are four parameters to consider within each size range. These are 
wavelength (X), size distribution parameter (a) (used in the distribution 
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law n(r) = Ar a ), and the two part index of refraction m = (n) - i (n 1 ). 
Thus each scattering phase function and associated cross section is a 
function in four dimensions P(0; a,X, n, n 1 ) and cr(a. A, n, n'). The empirical 
fit problem is thus to find a set of four dimensional functions to fit computed 
data to within some accuracy criterion. In addition the (0) dependence of 
(P) must be described for any chosen set of parameters. 

One approach would be to establish a four dimensional array of 
computed points and to interpolate between them. It was felt however that 
the search times involved in this procedure would be longer than desirable, 
and that the four dimensional interpolation procedure would be excessively 
complex. . 

Another approach, and the one taken here, is to fit, in a least squares 
sense, the computed data with continuous functions of the four variables. 
This has been done by applying a multivariate regression analysis to the 
computed Mie data for the extinction, scattering, and absorption cross 
sections (a^, cr g and a & ), and to the computed Mie phase function data at 
six selected angles, viz., 0°, 20°, 60°, 110°, 164°, 180°. The Mie data for 
each function (P,ct) was obtained for 81 points corresponding to all possible 
combinations of the three values assigned to each parameter, viz., a = 2, 
4, 7; X = 4000, 6000, 8000A; n = 1.3, 1.4, 1.5; and n' = 0, 0.02, 0.04. The 
continuous functions thus obtained for the phase functions and cross sections 
are listed in Appendix 5.2. 

There are several options for the choice of a "goodness of fit" 
criterion. In general one can force fit with a least path moment weight 
called the best Lp approximation (Kahng, 1972). When p = <» the result is 
the Chebyshev "maximum residual" fit. When p = 2 we have the usual 
least squares fit. For this problem some (p) in between would perhaps be 
desirable; however since a multivariate computer code was readily available 
only for p = 2, this fit criterion was used and the resulting functions are 
least squares fits in four dimensions. 
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Phase Function Interpolation 


To obtain P(0; a. A, n, n',) for (0) values other than 0^= 0°, 20°, 60°, 
110°, 164°, 180°, an interpolation formula is required. A least squares 
polynomial fit could be made to the selected angle points, however, 
comparable accuracy was found to be obtainable by fitting log-linear straight 
lines between the function values P^ at the selected angles 0^. Thus between 
the selected 0^, values of P(0) are given by 


p ( 0 il e i 6 i+i) “ P(6 i )exp 


In P (6 i+1 ) - In p (0^ 


( e i + i - e i) 


<e - e i> (3.3-1) 


Normalization 


The interpolation formula has the advantage that it can be easily 
normalized to insure proper conservation of scattering by shifting each 
modeled (P) up or down by an amount determined by a closed form integral. 
The constant of normalization is just 


5 {exp(m,0. ,)y. -exp (m. 6 . ) B . } 

c -= V 2 *.I P(0 i )exp(-m i e i ) ^ (3.3-2) 


i=l 


1 + m. 


i ‘ 


where 


y. = m iS in0. +1 -cos0 . +1 


(3.3-3) 


B . = m. sin0 . - cos0 . 
i i i i 


(3.3-4) 


and 


m. = (In P(0. +1 ) - In P (9 ± ) }/ (0 i+1 - 0 i) 


(3.3-5) 
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3.3.3 Error Sources 


It is possible to have considerable error in the phase function models 
for a particular set of parameters and still gain useful information from a 
model. The important thing is that the model represent an average over 
an ensemble of computed points that are representative of the phase function 
over the range of the parameters. For example if a modeled P(0) has an 
RMS deviation of 25% at some angle, but variations of a parameter (say 
n') over its range cause P(0) to change by say 200% then reliable partial 
derivatives and the concomitant inversion sensitivities can be obtained. 
The error sources and related limitations are described in the following 
subsections. 

Several strong trends in the data that should outweigh the RMS 
residuals are evident. These are noted below. 

1.0- IO/j Size Range 

1. There is a strong (n') dependence at both 164° and 180° as well as 
other angles in this region. The phase function is generally an order 
of magnitude lower for n' = 0.02 than for the n' = 0 case. The n' = 
0.04 case drops P(0) another 50%. 

2. The P(0°) values are twice as high for the n' = 0.02, 0.04 cases as 
for the n' = 0 case. 

3. (n 1 ) effects are also strong at 20° and 60° (200% to 500% variation), 
as well as at other intermediate angles. 

4. (X) effects are strong (300%) at 0° and 180°. As expected, the (X) 
dependence changes from P(0°) <* 1/X to P(180°) or x. 

5. (a) effects are as much as 400% at P(0°) and diminish with increasing 
(0)/ Beyond 110° complicated forms result. 

6. (n) effects are most noticeable and consistent at 164° (up to 300% 
variation). As (n) increases, P(0) increases. 


74 



0.1-1. Op Size Range 


1. The(n') dependence at 164° and 180° is strong; (200% to 500%) unless 
Q = 7, elsewhere it is very weak (10%). 

2. (n) effects are up to 200% at 164° and 180°, but diminish with increasing 
a and X . 

3. (X) dependence is strongly linked to (a) over all (0). 

4. (n) effects at 164° are strong (200% to 300%) unless a= 7, X= 8000A. 

5. (a) effects are strong at 0 (up to an order or magnitude), but correlated 

with (X), P(0°) °c 1 /c* A. " . 

0. 01.0. Ip Size Range 

1. (n 1 ) dependence is negligibly small over the range of all the other 
parameters. 

2. (n) dependence is at most a 10% effect at the larger angles for a = 2. 

3. (a) and (X) both strongly influence the phase function, but the effects 

are on a smaller (200%) scale than for the larger size ranges. 

Mie Code Errors 

To normalize the phase function generated by the weighted size 
distribution average, the Mie code performs a trapezoidal integration 
between the selected defining points, and adjusts the amplitudes by the proper 
constant to achieve conservative scattering. However, since much of the 
accumulated value of the normalizing integral 2 n / p(0) sin 0 d0 is obtained 
at small angles, the adjustment of the entire functional amplitude is very 
sensitive to the fineness of the angular defining intervals. This is 
particularly important for the larger particles which have rapidly decreasing 
phase function values near 0=0. In the (1.0-10p) size range, a decrease 
in the small angle (0) step size from 10° to 0.2° results in a factor of 
twenty adjustment of the phase function upward. The upwards adjustment 
is caused by the excessive integral accumulation under the trapezoidal 
function between 0° and the first defined point at 10° in the coarse sample. 
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Figure 3.3-1 shows an extreme case of this with angle step increments 
of 0.2° and 10°. Both curves in Figure 3 . 3-1 are normalized by integration to 
unity. The dashed curve is the actual curve under which the trapezoidal 
integration is performed rather than the computer- drawn straight line 
connecting the 0° and 10° points. In general the error due to the overall 
phase function adjustment for normalization is proportional to the difference 
between the scattering cross section as predicted, by the evaluation of the 
9=0° Mie series and the value obtained by integration of the phase function. 
For the (1.0-1 0>u ) size range the maximum value of this error is about 
25% forthe Ae intervals that were used in the calculations. In the (0. 1-1.0//) 
range the error is about 5% or less, and for the (0.01-0. Ip) range the error 
is less than 0.3%. 

Size Distribution Sampling Errors 

~Gt 

The power law distribution n(r) = Ar is sampled at discrete (r) 
values within the chosen size ranges. Each (r) choice corresponds to a 
particular Mie scattering result for a particle of radius (r), and the resulting 
phase functions are averaged by combining the individual (r) phase values 
with the (-Q!) weight. Since this process has the effect of averaging out the 
oscillations related to the Bessel and Legendre series within the Mie series, 
the results depend somewhat on the fineness of the Ar integration steps. 
Table 3.3-1 shows an example for the size range (1.0- lOp) with complex 
index of refraction. Differences between the phase function values for the 
range of 10 to 200 radius sample points are generally of the order of a 
few percent except for the 180° values which differ by 20%. In Table 3.3-II 
similar values are shown for non- absorptive aerosols (n 1 = 0). Thirty 
percent differences in the phase function values occur here fora few angles, 
and indicate that the number of radius steps required for a few percent 
accuracy is nearer 100 than 10. However the required computer time is 
directly proportional to the number of integration steps. 
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NO. OF STEPS IN INTERVAL NO. OF STEPS IN INTERVAL 


TABLE 3.3-1 
PHASE FUNCTION VALUES 



NOMINAL CASE, 

a=2, A= . 4 , n= 

II 

C 

• 

i — i 

n 

.02, 1-lOia 





PHASE ANGLE 





0° 

20° 

CD 

O 

0 

: O 

o 

164° 

180° 

10 

372 

.0367 

,00406 

.00120 

.00126 

.00122 

50 

368 

.0373 

.00411 

.00122 

.00118 

.00164 

100 

368 

.0371 

.00418 

.00119 

.00124 

.00152 

200 

368 

.0372 

.00418 

.00120 

.00123 

.00154 





TABLE 3.3- 

-II 





PHASE FUNCTION 

VALUES 





a=2, A= 

. 4 , n=l . 4 , 

n 1 =0 , l-10y 






PHASE ANGLE 




0° 

20° 

60° 

110° 

164° 

180° 

10 

256 

.205 

.0151 

.00264 

. 0114 

.0578 

50 

249 

.236 

.0190 

.00285 

.0134 

.0676 

100 

253 

.203 

.0196 

.00245 

.0141 

.0728 

200 

253 

.204 

.0196 

.00276 

.0138 

.0731 
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The derived aerosol model is based upon a 10 step integration, thus 
30% errors may be expected in the (1.0- 10a/) size range from this error 
source. Figure 3.3-2 shows a plot of the differences between the phase 
function values resulting from 10 step and 200 step size distribution 
sampling. The differences are expressed, as a percentage of the phase 
function value for the 200 sample case. It shows that the effects of changing 
the sample step size are essentially random with regard to angular position. 

Errors Related to Piecewise Function Fitting 

The errors just discussed were related, to the computational 
uncertainty of the phase function at a particular angle. The modeling 
procedure we have chosen connects the points f(P(0^), 0^) with exponential 
curves of the form P(0) = Aexp(B0); thus even if the parameter model P(0p 
n, n'.ar.X) were perfectly accurate there would be errors caused by the 
fluctuations of P(0) between the select model points P(0^). This would be 
true, whether the points were connected by exponential curves or by, e.g. a 
least squares polynomial fit. This error source would decrease in magnitude 
if the P(0^) sample points were. more closely spaced, however the six chosen 
0 ^ represent a compromise between accuracy and the need to limit model 
complexity. Each new 0^ requires a new Lp fit to find the modeled P(0p 
n, n\ a, X), however the data shows that the 0^ choices should be optimized 
for each size range separately. For example, Figure 3.3-3 shows the 
log-linear fit to the (0. 01-0 . lp) case havingthe most curvature . The fit is in 
error by a maximum of 14%. The RMS error for this case is 7%. For 
some of the better fits in this size range the RMS error is about 4%. The 
case illustrated in Figure 3 ,3-3 could be, improved by the addition of points 
defined around 30°, 90°, and 110°. However the choice of defining angles 
for a given size range is not a simple matter to be solved by examination 
of a single curve. The noticeable structure in the curves in the backscatter 
region moves over a range of (0) when the various parameters are changed. 

Figure 3.3-4 shows the six point fit to one of the worst-fit cases in 
the (0.1-1.0/y) range. In this size range the points of fit correspond quite 
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closely with the major points where curvature changes, however it can be 
seen from the figure that additional defining points would be helpful here 
at 130° and 176°. The maximum error for the (0. 1-1.0*/) case is 56%, 
with an RMS error of 21.5%, l<r. 

In the (1.0 -'1 0/y ) size range the standard points of definition do not fit 
the Mie curves well, and three or four additional points should be used for 
a better fit. A particularly bad-fitting case is shown in Figure 3.3-5. 
The minimum for the phase function is at 130° and is missed by the 
approximate fit point at 110° causing large errors at 130°. The 6-point 
function also ignores the secondary peak at 150°, and fits badly at 10°. 
Over the range of parameters considered, the position of the phase function 
minimum for the (1.0-1 0 a/ ) size range varies over the interval 110°-130°. 
Thus, improved modeling of this size range would require inclusion of extra 
defining points within this interval. The maximum error for this case is 
630% at 130°. The RMS error is 400%. 

The RMS errors listed in this section are based upon the sequence 
of discrete angles used in the Mie computations. Since these angles are 
defined with smaller increments at small and large angles, there is a built-in 
weighting which may adjust the RMS error value somewhat either upward 
or downward (depending on the particular fit) from the RMS error which 
would result from a uniform disbursement of sample points. 

In summary, the errors due to the piecewise log- linear fit at six 
preselected angles are moderate but within useable limits for the (0.01-0.1/u) 
and (0.1-1.0/u) size ranges. The errors are too large for the large particle 
range (1.0-10/u) which implies that an optimization procedure should be 
constructed for choosing the angular points of definition. 

Parametric Fluctuations 

It is of interest to examine the effects of small changes in the defining 
parameters (a.'X.n, n 1 ) to see whether or not the phase function is undergoing 
rapid fluctuations near the selected nominal parameter values. If this were 
the case there would be no guarantee that the value of P(6) at the nominal 
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parameter value was actually representative of P(0) in the region around 
the selected value. 

This problem was investigated by means of a set of Mie calculations 
with small perturbations from the standard points. Deviations from the 
case a = 2, X = 4000A, n = 1.4, n' = 0.02 were considered by running cases 
with a= 2.1, X = 4100A, n = 1.41 and n 1 = 0.021. The (1.0-10//) size range 
was considered because it exhibits the most erratic behavior and should 
represent a worst case. The results are tabulated in Table 3.3-III. 


TABLE 3, 3 -I II 


LINEA.RITY 

PARAMETER 

5 FOR A 

1 -lOy 

CASE 

Phase 

Angle 

(a) 

(n) 

(n') 

(X) 

0 ° 

5=2.3 

2.4 

0.12 

4.2 

20 ° 

9.0 

1.6 

0.2 

7.7 

60° 

0.36 

5.0 

2.7 

18.0 

110 ° 

0.91 

0.18 

0.27 

1.0 

164° 

1.2 

1.6 

0.12 

2.1 

180° 

0.39 

0.63 

0.027 

6.7 


. The number in the table is F = -7=— whe re AP and Ax aretheranges 

of the phase function and theparameter expressed as a percent of thenominal 
value. The quantities 6P and 6x are the perturbed phase function and 
parameter values also expressed as a percent or fraction of the nominal 
values. If all the functional behavior were linear these values would be 
unity. Typical changes in the phase function indicate the £= 4 corresponds 


80 



to quadratic dependence. The results show that at X = 4000A one has a 
rapidly varying phase function at 60°, however the excursion from the value 
at 4000A to the value at 41 00 A is a variation of 21% which is comparable 
to the scale of other errors for this size range. Otherwise the fluctuations 
are relatively calm. The small values for (n') are related to the fact that 
P changes rapidly in the forward and backward directions for very small 
deviations of (n 1 ) from zero. Thus the range of P is large between values 
corresponding to n' = 0 and n' = 0.02, hence the small values of (£) for 
(n 1 ). Since the modeling functions used were primarily quadratic, one expects 
to find relatively smooth behavior of the phase functions, with respect to 
parametric variations over the parameter range. 


Summary of Phase Function Modeling Errors 


The models listed in Appendix 5.2 exhibit a complex error structure 
in general. However one can estimate the maximum RMS errors by assuming 
that the various sources of error are independent and thus that the errors 
can be added. The error sources are, in approximately decreasing order 
of importance, regression residuals, piecewise fitting errors, , size 
distribution sampling errors, Mie Code normalizing errors, and higher 
order parametric fluxuations. 

The regression residuals are listed in Table 3.3-IV with both maximum 
and RMS values given. 

At the defined angles the total error consists of the regression error, 
the size distribution sampling error, and the normalizing error. If we 
take the typical size distribution sampling errors of Table3.3-I as la values, 
and the normalizing errors discussed above and combine them by means 
of a root mean sum of squares, the errors in Table 3.3-V are obtained at 
the defined phase function angles. 
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TABLE 3. 3. -IV 


RMS AND MAXIMUM RESIDUALS FOR PHASE FUNCTION PARAMETER MODEL 


Phase 

0.01 

-O.ly 

0.1- 

l.Oy 

1.0 

-lOy 

Angle 

o (%) 

Max. (%) 

o (%) 

Max. (%) 

o(%) 

Max. (%) 

0° 

0.3 

0.9 

3 

8 

3 

10 

20° 

o 

• 

to 

0.5 

1.8 

5 

11 

45 

o 

o 

VO 

o.l 

0.3 

4 

9.5 

8 

23 

110° 

0.4 

2.7 

4.5 

17 

18.9 

97 

164° 

0.25 

0.6 

6 

15 

31 

120 

180° 

0.2 

0.5 

10 

25 

32 

127 


TABLE 3.3.-V 

RMS ERRORS AT DEFINED ANGLES 
INCLUDING NORMALIZATION AND SAMPLING INTERVAL ERRORS 


Phase 

Angle 

0.01- O.ly 

0 .1— 1.0/u 

1.0-10y 

(%) 

(%) 

(%) 

0° 

2% 

6% 

25% 

20° 

2% 

5% 

27% 

0 

o 

VO 

1% 

9% 

35% 

110° 

3% 

15% 

32% 

164° 

4% 

9% 

43% 

180° 

4% 

16% 

46% 
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In between these "defined" values the piecewise fitting errors must 
be added in proportion to the distance from the defined angle. These errors 
have been discussed in one of the previous paragraphs, and since no attempt 
has been made to assess them, statistically they will not be combined with 
the above table to produce a grand model error table. For the (0.01-0. 1/u) 
size range the "in-between" errors are of the order of 10% at maximum. 
For the (0.1 -1.0a/) range the maximum fitting error is roughly 50%, and 
for the (1.0-10/y) range it is approximately 100%. A graphical display of 
the error sources is shown in Figure 3.3-6. 

3.3.4 Models of the Extinction, Scattering and Absorption Cross Sections 

The cross sections were modeled by applying the regression 
procedure described above for the phase functions. Since the cross sections 
for absorption and scattering sum to the extinction cross section, only two 
of the three cross sections need be modeled in each case. Thus the two 
with the smallest variances were chosen from each set of three models 
and used to produce the third by algebraic combination. The errors for 
these quantities are listed in Table 3. 3- VI. 

TABLE 3. 3 -VI 
CROSS-SECTION ERRORS 

0 .01-0 .ly O'. 1-1 . Oy 1.0-lQp 

a(%) Max ( % ) a(%) Max(%) a(%) Max(%) 

Error Error Error 

a (absorption) 3 10 25 29 

cl 

a g (scattering) 0.3 1 - - - 

a (extinction) - - 10 34 18 52 

e 

The standard errors for the unlisted cross sections can be determined 
if the correlations are known, however, since these errors were not used 
explicitly in the inversion procedure the correlations were not computed. 
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PHASE FUNCTION 



Fig. 3.3-1 Dependence of Normalization Adjustment 
Upon A 0 Sample Size. 
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Fig. 3.3-2 Percent errors due to size distribution sample intervals 
using 10 and 200 steps, (l-10y) aerosol size range 













PHASE FUNCTION 


Fig. 3.3-5 


Worst Case Fitting Errors, 

1.0-10y Size Range 
88 





0 . 01 - 0 . 1 ^ 



Fig. 3.3-6 Phase function error sources 
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4.0 SUMMARY 


4.1 CONCLUSIONS — AEROSOL PHYSICAL PROPERTIES 

The vertical distribution of aerosols from approximately 10 km 
upwards to 1 00 km can be determined globally by a satellite horizon inversion 
technique. Additionally, information such as the number density, the size 
distribution, and the complex parts of the index of refraction can be obtained. 


Three particle size ranges (0.01-0. 1//, 0. 1-1.0//, 1.0-10//) were 
considered in the simulations, and the results show that the quality of the 
physical characteristic information is size range dependent, i.e., that the 
relative observability of the physical characteristics changes from one size 
range to the next. For the smallest size range (0.01-0.1//) the number 
density (/>) and the size distribution parameter (ct) are observable, but the 
indices of refraction are not. However this result rests partly upon the 
assumption that the initial uncertainties of the four modeled physical 
parameters are as follows: cr = ±300%, cr ' = ±1, a , = ±0.005, 

cr no = ±0.05. Reductions in a and have the effect of increasing the 
observability of (n 1 ). A simple covariance propagation case in this size 
range, using four kilometer horizon sample intervals and five wavelength 
channels, reduces the initial variances as follows: ct 0o (±300%)-kt 0 = ±150%, 


=±0 - 1 ' °n - 


no' 


V 


a n’o‘ 


For the medium sized particles (0. 1-1.0//), using the same initial 
covariances, the simulations show that the complex part of the index of 
refraction becomes generally observable. The peak sensitivity of (n 1 ) jumps 
from a value of «0.01 to«0.4 which is in the highly observable range. A 
simple covariance propagationcase, using four kilometer horizon sample 
intervals and three wavelengths showed initial covariance reductions as 
follows: ^(±300%) -*v fi = ±60%, a ao (±l)-*’T Q = ±0.2, 

cr n , o (±0.005Ha n , = ±0.001, tr no (±0.05)-+or n = ±0.05. 


91 



All four physical parameters are found to be observable for the large 
particle (1.0-1 Op) size range. However the observability of the index of 
refraction is strongly dependent upon the scattering angle, i.e. the sun 
azimuth and zenith angles. For a scattering angle in the forward region 
(«50°), (n) is found to be unobservable. In the side- scattering region of 
the phase function(n') is unobservable and the observabilities of (n) and (a) 
are considerably reduced. However in the backscattering directions all 
four aerosol parameters are observable. 

An estimate of the state errors resulting from errors in the models 
was made based upon the filter gain and computed horizon intensity 
deviations. The calculations show that the aerosol modeling errors induce 
physical parameter estimation errors that are generally smaller than the 
final variances resulting from a covariance matrix propagation. However 
the error analysis also shows that a deviation of the albedo of 0.5 from 
the assumed value can cause large uncertainties in the state vector elements. 
Thus the albedo should be included as an element of the state vector. 

4.2 RESEARCH EXTENSION 

The results obtained in study show how well one can determine the 
parameters of a specific aerosol model. For example the inchoate model 
developed for this research uses a single parameter, power law size 
distribution. If the real aerosol size distribution is of some other functional 
form, then the results of applying the chosen model in an inversion would 
be to obtain a best power law fit to the real function.- If the best power 
law fit to the real function can be closely related to the real distribution, 
then the results are useful for aerosol research. To this end, then, it is 
important to attempt to generate the most realistic aerosol models possible. 
A first step in this direction is to insure that the ranges of the chosen 
variables cover the real physical range. This has been done in the model 
used here in every respect except the aerosol size range. Hence one of 
the next steps in a continuation of aerosol model development would be to 
reconstruct the model with variable particle radius limits. This would 
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increase the size of the state vector by two elements, however, the present 
model could be used to optimize the selection of wavelength channels. By 
fixing wavelength, the dimension of the model could be reduced by one 
element for a net increase of one dimension in the model. 

To insure that fitting errors are minimized in extended aerosol models, 
a procedure should be constructed which determines a minimum set of 
phase function defining angles consistent with accuracy goals. The models 
should also make maximum use of the analytical relationships derived from 
ray optics (Van de Hulst, 1957) and information theory (Whitney, 1972). 

The design of the filter algorithm should be further extended to include 
the additional states resulting from improvements in the aerosol model 
and from the inclusion of polarization measurements. The effects of model 
type error sources such as approximation errors in the radiative transfer 
simulation, fitting errors in the aerosol models, and relative fluctuations 
in the solar power output should be explicitly included in the filter 
formulation. As the aerosol models become more complex, additional 
sources of information will be required besides the measurements at 
different wavelengths. Two possible sources are measurements at different 
solar angles and measurements of polarization. The present radiative 
transfer simulation for a curved atmosphere includes the necessary 
polarization calculations, however, these simulations would have to be 
checked against other results to insure their accuracy. Accurate inclusion 
of cloud scattering effects at low sun angles (<10°) would require 
modifications of the existing code as well as simulations of the transfer 
problem within a cloud. 
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5.0 APPENDICES 

5.1 THEORETICAL HORIZON PROFILES 

This appendix contains the theoretical multiple scattering horizon 
profiles discussed in Section 2.3. The profiles illustrating the effects of 
changes in various measurement parameters are arranged as follows. 

Figure Numbers Variable Measurement Parameter(s) 

5.1- 1 to 5.1-12 Season; for 12 combinations of latitude 

and wavelength. 

5.1- 13 to 5.1-28 Wavelength; for 16 combinations of cloud 

cover, ground cover and zenith angle. 

5.1- 29 to 5.1-64 Zenith angle; for 36 combinations of azimuth 

angle, wavelength, and ground albedo. 

5.1- 65 to 5.1-72 Ground albedo; for 8 combinations of 

zenith angle and wavelength. 

2 

The units of solar irradiance (SUN PWR) are in /uwatts/cm A. 

Figures 5.1-73 to 5.1-84 show the vertical distributions of neutral 
and ozone densities for 12 combinations of season and latitude. These 
figures illustrate the mean and one sigma standard deviation values about 
the mean. 

Neutral density data was obtained from Valley, 1965 and from the 
data compiled by Groves, 1970. The designations winter, spring, summer, 
and fall refer to the months January, April, July, and October. The oz;one 
density data was generated from the modified Fermi distribution function 
developed by Wu, 1970 and the peculiar appearance of the ozone density 
distribution for the minus one sigma curve is due to the fact that the one 
sigma values are equal to or greater than 100 per cent of the mean value 
below approximately 10 km and again above approximately 50 km. Aerosol 
extinction data was compiled from a variety of sources by Malchow, 1971 
and the composite model illustrated in Figure 5.1-85 was constructed. 
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Latitudinal variations in the aerosol extinction data were accounted for on 
the basis of data gathered by Salah, 19 71 as illustrated in Figure 5.1-86. 
Due to a lack of sufficient seasonal distribution, this aerosol extinction 
data was used to represent the four seasons. 
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-13 Variable Wavelength Horizon Profiles 

(Zenith=30°, Cloud Top at 4 km, Albedo=.4) 
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-14 Variable Wavelength Horizon Profiles 

(Zenith=30°, Cloud Top at 4 km, Albedo=.8) 



C.K.W. HORIZON PROFILES .. f DIFFERENTIAL HORIZON PROFILES 
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Fig. 5.1-15 Variable VJavelength Horizon Profiles 

(Zenith=30°, Cloud Top at 7 km, Albedo=.4) 
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Fig. 5.1-18 Variable Wavelength Horizon Profiles 

(Zenith=30°, Cloud Top at 10 km, Albedo=.8) 
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-21 Variable Wavelength Horizon Profiles 

(Zenith=80°, Cloud Top at 7 km, Albedo=.4) 






-23 Variable Wavelength Horizon Profile 
(Zenith=80°, Cloud Top at 10 km. A1 
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Fig. 5.1-25 Variable Wavelength Horizon Profiles 

(Zenith*=0°, Cloud Top at 0 km, Albedo=.3) 




oooos 



ro 


O 

(0 a) 

Q) XI 

r— ! rH 

m 

o ^ 

m e 

CM r* 


c 

o 

N 

-H 

M 

O 

a 

X! 

«P 

tn 

C 

0) 


•p 

fd 

a 

o 

Eh 

H3 

0 

O 

iH 

U 


a) 

> 

<d o 

£ o 


0) 

rH 

X 

<d 

-H 


00 

II 

•C 

■P 

•rH 

<u 


(d CS3 
> ~ 


VD 

CN 

I 


lD 


Cp 

•H 

P4 


122 





> 50 00 . 




Z 

V 


LlI 

o 

D 

H 

P 

< 


z 

v 

Ixl 

Q 

D 

t 

P 

< 


00 

li 

O 

Tl 

W cu 

QJ XJ 

iH r-H 

UH 

0 * 

u g 
cu ^ 


G 

o 

N 

u 

o 

ffi 

-P 

tr» 

G 

<U 


-P 

CU 

o 

E-J 

G 

O 

rH 

u 


d) 

> 
ftJ o 
IS o 


Q) 

i— I 

\Q 

ftf 

•H 

P 


00 

II 

& 

-p 

•H 

G 

<D 


fd [S3 
> ~ 


00 

(N 

I 


LO 


to 

'rH 

Cm 


(ans v ? wd /siivnrt) 30Nviava 

0 2 


124 




o 





(x) 30NVH3 A11SN31NI 


S S S o 8 

I 



o o o o o o 

f-4 rH T“^ r-4 r*4 


(U31S V ? WD /Slivnrt) 


WINTER 

LATITUDE = -*0 DEO. NORTH 



(Azimuth=0 , A=3000 A, Albedo=.2) 






































2nith 








o s s 

I 





u 

< 


*H 


m 


o 

in 

n 

• 

Pu 

II 


0 

G 

vd 

0 

0) 

N 

X 

•H 

1 — 1 

H 

< 

o 


ffi 



cc 

a) 


iH 

o 

CP 

o 


% s 


XJ r< 

c o 

<1) O 

tsa cn 
II 

<D jC 
rH -P 
XI 3 

* e 

•H -H 
U N 
rd i< 
> w 



CM 

I 


in 


tp 

•H 

P4 


(U31S V 7 WO /SllVttrt) 3DNViava 

0 « 


138 







DIFFERENTIAL HORIZON PROFILES 
VARIABLE ' ZENITH ANGLE (DEGREES) 
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SPRING 40N LAT f t SPRING ^ON L£T 



Fig. 5.1-78 Neutral and Ozone Density Profiles (Spring 40 N. LATITUDE) 
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Fig. 5.1-80 Neutral and Ozone Density Profiles (Fall, 40 N. LATITUDE) 
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Fig. 5.1-83 Neutral and Ozone Density Profiles (Summer, 70°N. LATITUDE) 
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Figure 5.1-85 Aerosol Extinction Profile 
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5.2 MODELED AEROSOL FUNCTIONS 

This appendix gives the analytic expressions obtained (by a 
multivariate regression procedure discussed in Section 3.3 for scattering 
phase-functions P((h; Q, X, n, n') at A = 0°, 20°, 60°, .110°, 164° and .180°, 
the total cross-section cr^a, X, n, n') and either the absorption cross-section 
a (a, X, n, n') or the scattering cross-section a (a, X, n, n'). These 

o. S 

expressions are grouped according to the particle size ranges (0.01-0.1//, 

0. 1-1.0//, and 1.0-10//). 

In the 0. 01 to 0. 1// range we have: 

P(0°) = exp(-2.937E-2 + 5.234E-2 aX - 4.747 X - 5.812E-2 aX 2 

- 3.388E-2 a 2 + 3.259 X 2 + 6.963E-2 a 2 X + 1.195E-1 n 2 

- 3.961E-2 a 2 X 2 - 3.744E-1 nX) 

r\ 

P(20°)= exp (-2 . 396E-1 + 4.447E-2 aX - 4.687, X - 4.803E-2 aX 

- 3.229E-2 a 2 + 3.069 X 2 + 6.728E-2 a 2 X + 1.295E-1 n 2 

- 3.887E-2 a 2 X 2 - 1.451E-1 n 2 X) 

P (60°) = exp (-1.776 + 6.964E-2 aX - 1.906 X - 9.225E-2 aX 2 

-1.402E-2 a 2 + 1.388 X 2 + 2.051E-2 a 2 X + 9.914E-2 n 2 

- 2.343E-1 nX - 4.050E-3 na 2 X 2 ) . 

P(110°) =exp(-4.482 + 4.928 X + 3.946E-2 a 2 - 3.747 X 2 - 1.013E-1 
a 2 X - 2.372E-1 n 2 + 6.485E-2 a 2 X 2 + 2.946E-1 nX - . 

+1.534 E-l n 2 X + 1.100E-2 n 2 X 2 a) 
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P (164°) =exp (-6.612 - 1.685 aX + 1.227E+1 X + 1.166 aX 2 

- 8.577 X 2 - 3.798E-1 n 2 - 7.140E-3 a 2 X 2 + 4.126E-1 n 2 X 
+ 5.990E-1 a + 5.850E-3 na 2 ) 

P (180°) =exp (-6.686 -1.711 aX + 1.249E+1 X + 1.179 aX 2 

'-'8.693 X 2 - 3.715E-1 n 2 - 6.800E-3. a 2 X 2 + 4.027E-1 n 2 X 
+ 6.117E-1 a + 5.800E-3 na 2 ) 

a = exp (1 . 511E+1 - 1.779E-1 aX +1.189E-2 naX 2 +3.959 X 2 
s 

- 1.604 a - 9.541 n” 1 - 1.043E+1 X + 2.689E-2 a 2 ) 

a = .007259 + 77.895 n'{ (2.58 (a/2) 2,35 - . 79a) (X/. 4) 1 * 3 } -1 
a 585 

(1 + .0208 X(a-2) (7-a) { | (n-1.234) | / . 066 }“ (X/ * 4) * 

- {.1667 (a-2) (a-4) (n/1 . 3) 5 ’ 75 ( . 06-n ' ) H 1+7 . 5 (X- . 4) (X— .6) }} 
In the 0.1 to ly range we have: 

P(0°) = exp (6 . 415E-1 + 2.302 naX + 1.992E-1 nX 2 a + 5.082 nX 2 

+ 1.841E-1 a 2 X 2 + 9.519E-1 na - 1.404E-1 na 2 + 9.720 nn' 

- 1.533 nn'aX - 1.560E+1 nX + 1.679E-1 (na) 2 

+ 2.865E-1 nn'(aX) 2 - 2.336 aX 2 + 3.003 n 2 - 1.569 n 2 a 
+ 9.547. X - 2.570E-1 na 2 X + 3.095E-1 n 2 aX - 6.137 nn'X 2 

- 9.415E-1 nn') 

P ( 2 0° ) = exp ( - 3 . 6 6 8 - 9.038E-3 (na) 2 + 2.348E-1 naX - 3.470E-3 

(na) 2 + 2.986 n" 1 ' - 8.749 a 2 X + 4.534E-2 (naX) 2 - 4.293 X 2 
+ 1.473E-1 na - 3.111E-1 (nX) 2 a + 7.388 nX 2 + 1.166E+1 
Xn'n 2 - 2.672E+1 n'X - 1.666E-1 n'(Xna) 2 - 2.416 (nX) 2 
+ 2.292E-1 n'(Xa) 2 + 8.322E-2 a + 2.412E+1 X 2 n' + 1.235E-1 
n'n 2 a - 9 . 250 (Xn) 2 n ' ) 
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P (60°) = exp (-5 . 344 - 1.100E-2 n 2 a + 3.714 n 3 X + 1.168E-1 n (aA) 

0 o 2 

- 2.729E-2 (naX) ^ - 4.697 n X + 2.213E-1 (n) n a + 2.068E+1 

n'X/a - 1.715 X n + 2.363 naX - 1.429 n 2 aX - 4.949E-1 n n 
+ 3.470 n 2 /{ (l+5n ' ) a} - 2.016E+1 n 2 a -3 + 5.328E-1 (n')‘ 5 
+ 9.547 n 2 a“ 2 - 3.757E-1 n 2 (l+5n') -2 - 2.719 n' + 1.058E-1 
exp (-10 n'a) - 4.793E-3 (na) 2 X - 2.276E-1 n 2 + 2.996E+1 
(n'X/a) 2 + 4.266E-2 n 2 (l+5n')~ 1 ) 

P (110°) =exp (4.751 - 4.099 n 2 aX - 8.511 n + 3.866 aX + 9.015E-2 
(naX) 2 - 3.342E+1 n* + 8.520 Xn'a - 3.278E-12 na 12 /X 
+ 6.585E-1 a 2 X + 2.164 n 2 a - 5.439 aX 2 - 3.228E+1 nX 2 

- 9.783E-2 nn'(aX) 2 - 7.001E-1 (aX) 2 + 4.486E+1 nX 

+ 3.729E+1 X 2 - 1.196E-1 na 2 - 5.399E+1 X + 1.259E+1 nn' 

- 3.436 nn'aX 2 - 5.098 n'(na) 2 + 1.143E+2 X(n') 2 + 6.578 
naX 2 - 2.303 na - 2.545E+1 a(n'X) 2 ) 

P (164°) =exp (2 . 472E+1 + 7.151E-1 (a-2) (a-4) (X-. 4) (X/. 6) " 1 * 08 

+ 2.546 (n-1.35) |a-7| (n/1.4) _3,88 - 4.058E+1 nn' + 7.673E-1 
nn'(aX) 2 - 1.056 a 2 - 5.101E+1 X 2 - 1.843 a 2 X + 7.514E-1 
nn'a 2 X - 1.471 a + 1.769E+1 aX - 7.267E-1 (1+50 n')/a 
+ 2.921 na 2 + 5.524E-1 X(na) 2 + 7.434 (naX) 2 - 1.854 
n (aX) 2 - 2.226 X/(l+50 n') + 4 . 555E+2 ,n (n ' ) 2 + 3.278E+1 
nX 2 - 7.245 aXn 2 - 1.336E+1 na - 1.815E+1 n + 9.171 n 2 a 

- 1.379 (na) 2 + 1.352 (aX) 2 + 2.562E-1 a(nX) 2 ) 

2 

P (180°) =exp (1.340 + 5.193E-2 ln{ (a-2/X) (2/n ] (1+50 n’) -1/X } 

- 2.180E+1 n - 1.070 nn'(aX) 2 - 1.092E-1 (na) 2 - 1.379E-1 
na 2 + 3.134E-1 Xa 2 + 3.732 nn'aX 2 + 5.469E-1 (naX) 2 
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- 3.927E+1 nn' + 8.804 n'aX - 8.464 n 2 aX+ 1.187E+1 X 2 

- 3.627E+1 nX 2 + 3.261 aX - 1.255 (aX) 2 + 1.495E-1 
n'(na) 2 + 1.427 n' + 3.153 naX 2 + 3.236E-1 (a-2/X) 2 
+ 4.157E+1 n (n ' ) 2 + 3.377 aX 2 + 8.490 nX + 3.580 an 2 

- 5.849E+1 X) 

4 -i 4 • -1 

a_ = exp{2.3{X A.a^ - exp ( B . cP ) } { 1+ . 036a (n-1 . 4) } -1 } 

1 j=0 3 j=0 3 

where A Q = -8.896 E-l A^ = 1.288 

A 2 = -5.316 E-l A 3 = 6.727 E-2 

A. = -2.855 E-3 
4 

B q = 1.948 B^^ = 1.159 E-l 

B 2 = -2.551 E-2 B 3 = 2.491 E-3 

B 4 = -8.889 E-5 

a = -1.518E-3 + FX {1 +12.5 (X-.6) (X-.8) (.215 + DX + .-5 (n-1. 3) 
a 

CX + (n-1. 3) (n-1. 4) BX> + 25 (X- . 4) (X- . 8) ( . 25- . 03 (7-a) 

- .7 (n-1. 3) ) - 12.5(X-.4) (X-.6) (.47 - .032 (7-ct) - .7(n-1.3) 
EX - .2(n-1.3) (a-4) (a-2)/15) } (1 .61 (a/2) 2 ‘ 32 - .61) _1 

where FX = 1.372 (n'/.02) AX 

AX = .585 + 2.5 U-.4) (.138 - .0163(a-2)) + .0644 (a-2) 
CX = | (a-2) |a“ a//8 
BX = a 2 CX/8 

DX = | (n-1. 3) | (n/1.4)" 10 * 21 
. EX = | (a-7) | a/12 
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In the 1 to 10y range we have: 

P(0°) = exp (8.457 - 1.357 a + 6.576 X + 9.534E+1 n' - 6.383 a 2 

- 2.374E+4 (n') 2 - 8.454 X 2 - 1.596E+1 Xn'a 2 + 2.076E+3 
X (n ' a) 2 - 7.498E-1 aX 2 - 1.161E+1 (aX) 2 + 9.258 na 2 

+ 8.853 X (na) 2 - 3.303 (na) 2 - 3.573E+1 nn' + 1.211E+2 
Xn'a + 1.375 Xn'n 2 - 6.465E+1 Xa(n') 2 - 3.425E+3 Xn(n') 2 

- 9.796E+1 (Xan') 2 - 3.656E+1 an' + 1.262 n 2 + 5.706E+1 

na(Xn') 2 + 4.348 n'a 2 + 5.356E+2 a(Xn') 2 + 2.365E+3 
n(Xn') 2 + 1.036E+3 X(nan') 2 - 7.508 n'(Xna) 2 - 6.004 
(naX) 2 + 1.221 aX + 1.677E+1 n(aX) 2 - 2.909E+3 nX(an') 2 

- 1.155E+4 (nn') 2 + 3.359E+4 n(n') 2 - 1.199E+2 an'X 2 

+ 9.740 nX 2 + 1.736E-3 nn'a 2 - 1.195E+1 nX - 2.378E+1 
nXa 2 + 1.717E+1 Xa 2 + 2.148 nn'(aX) 2 ) 

P(20°)= exp (-6.551 - 5.221E+1 nn' - 6.469E+2 naX(n') 2 + 6.866E+1 
n(aXn') 2 + 3.327E+2 nn'(aX) 2 - 8.600E+4 (nn') 2 
+ 3.156E+3 aXn ' - 1.193E+2 n'(anX) 2 - 2.377E+2 n'(aX) 2 
+ 7.901/n + 2.525 n 2 - 5.764E+1 n' - 1.585E+3 X 2 

+ 5.231E+2 Xn 2 + 9.972E-2 (na) 2 - 1.822E+1 Xa 2 

+ 2.823E+1 (aX) 2 - 8.998 X(na) 2 + 2.520E+1 nXa 2 

+ 2.304E+3 nX 2 - 8.260E+2 (na) 2 - 1.941E-2 na 2 + 2.430E+5 

n(n') 2 - 1.689E+5 (n') 2 - 3.983E+1 n'X 2 - 1.402E+1 nan' 

+ 9.910 a (nX) 2 + 1.379E+1 (naX) 2 - 4.360 naX+ 1.000E+3 X 

- 3.913E+1 n (aX) 2 - 1.560 na + 1.587E+3 Xan'n 2 - 4.434E+3 

nn'aX - 1.463E+3 nX - 1.903E+1 naX 2 + 1.422E+1 aX + 6.747 

an'n ) 
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P (60°) = exp (-2 . 065E+2 - 1.942E+3 n' + 1.348E-1 aXn 2 + 3.908 

(aX) 2 + 5.889E+2 an'X 2 - 2.647E+1 (n') 2 + 3.968E+2 a(Xn') 2 
+ 1 , 789E+2 n -1 + 3.761E+1 n 2 + 3.453 aXn' + 2.123E+2 nX 

- 3.399E+2 na(n') 2 - 2.754E+1 an 2 - 1.225E-2 a 2 - 8.192E+1 
Xn 2 - 1.354E+2 X + 5.389 Xa 2 - 4.490E+2 na(Xn') 2 

- 1.054E+3 n'n 2 + 7.534E+1 na + 1.800 (naX) 2 + 1.274E+3 
(nn') 2 + 3.513E+2 an'(nX) 2 - 2.453 n'(Xna) 2 + 3.865E+1 
n(aXn') 2 + 1.520 an'n 2 - 4.816 aX 2 - 5.124E+1 a - 5.320 
n(ctX) 2 + 2.786 E+3 nn' - 9.080E+2 nn'aX 2 - 5.976E+2 

X (nn ' ) 2 + 1.106E+1 (ann') 2 - 7.938 nXa 2 + 2.916 X (na) 2 
+ 3.193 naX 2 ) 

P (110°) =exp (2 . 558E+2 + 5.608E+1 a(nX) 2 + 1.412E+4 n' - 2.236E+4 nn' 
+ 8.719E+3 n'n 2 + 9.753E+3 (n') 2 - 2.447E+2/n - 4.403E+1 n 2 

- 1.587E+2 naX 2 + 2.201 a(nn') 2 - 1.782E+1 aX - 1.075 an 2 

- 1.409E+4 n'X 2 + 1.440E+4 X(nn!) 2 +.1.U3E+2 ctX 2 

- 4.937E+3 Xn'n 2 - 1.838E+1 an'n 2 - 2.539E+2 aXn' + 7.063 

nn' X + 1.474E+4 nn'X 2 - 3.422E+3 n' (nX) 2 + 1.419 (ann') 2 
-5.161 Xn' a 2 - 5.923 n' (na) 2 4.134E-3 (na) 2 - 2.600E+3 

an'X 2 - 5.367E-1 X(ann') 2 - 3.557E+3 (nn') 2 - 2.589E+4 
nX(n') 2 - 7.125 aXn 2 + 2.336E+1 naX + 6.357E-3 a + 1.431 
na + 5.565E+3 (Xn') 2 - 1.618E+3 an'.(nX) 2 +.4.130E+3 ann'X 2 

.+ 1.549E+2 aXn'n 2 + 1.079 nn’a?) 

P (164°) =exp (2 . 689E+2 + 1.618E+4 n' + 8.978E-2 a(nX) 2 -2.602E+2/n 

- 2.900E+4 (n 1 ) 2 - 3.299E+3 an'X 2 - 3.436E+4 an(n'X) 2 

+ 6.153E+3 n'n 2 - 2.036E+4 nn' - 4.332 n'(aX) 2 + 1.493 aX 
+ 1.669E+4 (nn') 2 + 2.443E-2 na 2 + 2.274E+3 nn'aX 2 . 
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+ 4.619E+1 cm' (nA) 2 - 5.877E+1 na - 1.269 E-l (aX) 2 

- 1.429 An 2 - 9.151E+3 aA(nn') 2 + 7.659E+1 A(ari') 2 

- 1.297E+3 aAn'n 2 - 8.899E+4 n(An') 2 - 5.749E+3 An' 

- 4.313E+1 n 2 - 8.996E+3 a(nn') 2 + 2.629E+3 aAn' 

+ 1.798E+4 a (n' ) 2 - 1.885E+3 an'? + 1.339E+3 ann' 

- 1.059E+5 aX(n') 2 + 5.003E+4 a (An',) 2 + 8.643E+4 aA(n') 2 
+ 4.100E+1 a + 1.062E+1 an 2 + 4.159E+3 Ann' ■+•' 1.225E+5 
(An') 2 ) 

P (180°) =exp (1 . 759E+2 .+ 1.334E+2 n' + 3.142E+4 aA(nn') 2 - 2.056E+3 
A ( an ' ) 2 - 3.018E+1 n 2 - 1..074E+3 An'a 2 •+. 4.978E+3 (n') 2 
+ 5.162E+2 aAn'n 2 - 2.187E+1 (aA) 2 - 1.068E+1 (naA) 2 
, + 6.399E+2 rin'Aa 2 + 2.924E+5 naA(n') 2 + 4 .239E+1, An' (na) 2 

+ 1.371E+2 aX +6.876 aXn 2 - 9.354E+4 a(n') 2 + 5.932E+2 
• A(nn') 2 - 6.690E+2 nn' (aA) 2 - 1.508E+3 An' - 1.663E+2/n 

- 1.945E+2 naA - 6.090 n'an 2 + 1.004E+3 n' (aA) 2 - 9.508E+3 
a (nn 1 ) 2 + 8.024E+4 an(n') 2 + 1.169E+3 n'A 2 + 1.629E+5 

na (An' ) 2 - 2.267E+5 a(An') 2 + 3.438E+5 aA(n') 2 + 3.062E+1 
n (aA) 2 - 3.701E+2 an'(nA) 2 + 1.939E+2 n'a 2 - 1.748E+2 nn'a 

- 8.555E+1 n' (na) 2 + 1.618E+3 nA(an') 2 ) 

a T = exp (-1 . 498E+2 - 1.020E+1 a + 3.730 a 2 + 8.622E+2 nA 

- 2.967E+2 n ' (aA) 2 + 1.548 (na) 2 - 5.009 na 2 - 2.402E+2 

n' (naA) 2 - 4.761E+2 n'Aa 2 + 1.484E+2 (nA) 2 + 4.393E+1 n'a 2 
+ 3.614E+2 X 2 - 6.377E+2 A + 2.016E+1 n 2 - 4.702E+2 nX 2 

- 2.852 Xn 2 + 5.636E+2 nn' (aA) 2 + 1.086E+1 n'(na) 2 

- 3.906E+1 nn'a 2 + 1.545E+2/n - 1.476E+2 An'(na) 2 
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+ 5.258E+2 nn'Xa 2 + 8.306 na + 3.476E+3 an’ (nX) 2 

- 5.662E+2 aXn'n 2 - 1.181 aX 2 + 3.284 a(nX) 2 +. 1.239E+1 aX 

- 1.423E+1 naX + 5.827E+3 an'X 2 - 9.137E+3 nn'aX 2 
+ 9.489E+2 nn'aX - 8.068E+1 n'a + 3.763E+1 n'n 2 

- 1 . 981E+3 rXn ' +. 2.406E+3 nn' X - 1.159E-1 X (na) 2 

- 6.879E-1 (aX) 2 + 8.424E-1 nXa 2 - 7.663E+2 Xn'n 2 
+. 8.076E+1 (Xn' ) 2 ) 

a = 1.195E-3 - 2.507E-4 a + 1.087E+3 n' + 3.008E+3 a(n') 2 

+ 2.789E-5 a 2 - 2.547E+2 (ah ' ) 2 + 1 . 191E+1 n'a 2 - 1.091E+4 
(n') 2 - 1.940E+2 an' - 6.414E-4 Xn 2 + 4.157E+3 X(n') 2 

- 3.911E+2 Xn' - 6.612E+2 nn' + 1.380 nan' - 4.455E+1 
a (Xn' ) 2 + 2 . 947E-2 n' (aX) 2 - 2.262 n' (na.) 2 - 2.715E+2 

h' (nX) 2 + 2.815 an'X 2 + 1.853E+2 n'n 2 + 2.201E+1 X(an') 2 
+ 2.609E+2 nn'X + 5.240 nn'a 2 - 3.396E+2 n'X 2 - 2.056E+2 
aX(n') 2 ; - 5.251 nX(an') 2 + 1.317E+3 n(n') 2 -r* 1.977E+3 
nX(n') 2 
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5.3 MULTIPLE SCATTERING CODE C.K.W. 

5.3.1 Introduction 

Because it has not been given more than a cursory description in 
print, and because its speed makes the model concept potentially valuable 
to the scientific community as a whole, and because it is the model to be 
adapted for all our future work we present in this and the following 
subsections a brief description of the multiple scattering code C.K.W. in 
it's present state of development. 

For the program C.K.W. the atmosphere is modeled as a set of 
concentric spherical layers, each characterized by extinction coefficients 
for Rayleigh scatterers, aerosols, and ozone, and by an index of refraction, 
all ‘assumed constant within a given layer. The acceptability of assuming 
constancy depends on the thickness of layers. Typically, we use 50 to 60 
layers of 2 km thickness, and the calculated horizon profile is insensitive 
to further refinement. 

The calculations naturally fall into a hierarchy of increasing difficulty. 
The simplest problem is single scattering, which is handled in great detail, 
with the full spherical geometry, exact integrations including refraction, 
detailed phase function data with polarization, and a Feugelson law for 
albedo (Kondratyev, 1969 ). The next case includes some nearly forward 
scattering on the path from the sun to the site of one major scatter into 
the detector line of sight. This is done with almost as much refinement as 
single scattering. The next, more difficult case includes some nearly forward 
scattering on the path to the detector. For this case, an approximate 
estimation of coupling between scan lines is introduced. The final case is 
full multiple scattering, which requires some additional approximations 
for the geometry. The following sections detail the methods and 
approximations used to treat each of the successively more difficult cases. 
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5.3.2 Single Scattering 


Basic Geometry 

The modeling of the horizon profile requires that we first of all 
establish a coordinate frame in which all the vectors that occur in the 
calculation can be expressed unambiguously. For this purpose, we affix a 
coordinate frame to the detector. The (z) axis is along the local vertical 
at the detector, the (x) axis is perpendicular to the plane of scan, and the 
(y) axis lies in the scan plane, pointing roughly in the direction from which 
the detector seeks light (Figure 5.3-1). This choice of coordinate frame 
is arbitrary, but convenient on several counts. First, this coordinate frame 
is independent of the detector altitude. This is useful because one of the 
applications for horizon profile modeling is in spacecraft navigation, for 
which the altitude may be unknown. Second, the zenith and azimuth of the 
sun can be measured directly from the craft, so the vector representing 
the input sunlight can be defined directly for this frame. Third, in the 
special but common case that the sun lies in the scan plane, the (x) axis is 
perpendicular to the scattering plane, so there is a coincidence in regard 
to the naming of polarization states: both the usual geometric convention 
and the convention used in scattering theory agree that polarization along 
(x) corresponds to the same direction in Poincare space. 

A second geometric requirement to be met is the specification of 
the earth and its atmosphere. Thus we introduce the height of the detector, 
(HDET), and the earth radius (RE), shown in Figure 5.3-1. Also we fix 
the radii to the bottoms of each of the layers, R^: 

^ = RE + L KM (5.3-1) 

where (KM) is the single-layer thickness and (L) is an integer member of 
set (NL) defining the number of layers. Next we precalculate geometric 
quantities that will be used over and over in the horizon profile calculation. 
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These are the scan path increments through the layers, DX^ for L = 0(1)NL. 
For the scan path tangent to the earth, the path increment through layer 
(L) is 

• DX l £ /2 RE KM (/L+I - Vh) (5.3-2) 

Almost these same path increments occur on all the other scan paths too, 
and so the numbers are saved for use on all the scan paths. 

The final item of geometric information to be handled is the location 
of any cloud. We denote the altitude of a cloud by the variable SOLIDTOP 
in Figure 5.3-1. If the cloud happens to fall between layers, then it occludes 
aportionof some layer. This portion is called SOLIDTIP. In the subsequent 
calculations, any path increment through the partially occluded layer should 
be adjusted from the appropriate DX^ above. The adjustment factor is 


Average Attenuation 


SQ 


/KM - /SOLIDTIP 
/KM 


(5.3-3) 


The attenuation within a path increment (DD) is one of the factors 
that determine how much radiation is scattered to the detector from (DD). 
At any point in (DD), the radiation arriving from the sun is damped by 
passage through a total optical depth denoted by (TTS). Similarly, the 
radiation scattering to the detector is further attenuated by passage through 
a total optical depth denoted by (TTD). Thus at a point, the attenuation is 


-TTS - TTD 
e 
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A full contribution to the output from the whole increment (DD) is of the 
form 

/ e" TTS “ TTD ' Z a. P .(t|») dx 
DD i 

where a. and p.(0) are the scattering coefficient and phase function for the 
th ^ ^ 

i constituent, and (dx) is. a differential increment of length. Assuming 
constancy of all extinction coefficients in the layer (or at least their local 
ratios), we split the contribution into two factors: 

1 / -TTS - TTD , 

DD J e dx 

DD 

and 

X a i Pi^) DD 

The first factor is the average attenuation, and its evaluation is the subject 
of this subsection. 


To perform the integration in the average attenuation we introduce a 
change of variables: . 


dx = J d (TTS + TTD) 

Here (J) is the Jacobian of the transformation, 

J = 1/ (d (TTS + TTD) / dx) 


(5.3-4) 


(5.3-5) 
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We assume the derivative in (J) can be approximated by 

. ■ A (TTS + TTD)/DD 


Thus we reexpress the average attenuation by 


' 1 • 

A (TTS + TTD) 



-TTS 


TTD 


d (TTS 


.+. TTD) 


This can be evaluated immediately as 


DAMP = -Ae“ TTS TTD /A (TTS + TTD), , (5.3-6) 


The above trivial result is modified slightly by the inclusion of the 
phenomenon of refraction. The basic principle involved is that light travels 
not exactly in a straight line, but rather in a slightly curved line. At a 
point on the line, we let (n) denote the refractive index, (r) the distance to 
, the earth center, and sine the angle of incidence on a surface parallel to 
the earth. The curved line is described by 4 


(n)(r) sin0 = constant 


(5.3-7) 


where the constant is the nominal tangent height Rjj of the scan line. 
Referring to the greatly exaggerated Figure 5.3-2, we see that where the 
unrefracted path would be intersecting the bottom of layer (L) at an angle 
whose sine is R^/R^, the refracted pith is passing through a somewhat 
lower point and at a lesser angle. Let (y) be the nominal scan line, and 
(z) be the perpendicular displaceinent of the true path from the unrefracted 
one'. Then the true angle is 


0 -y arcs in 

Of 



arcsin 


^dz 

& *L 


(5.3-8) 
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We also have 


6 y arcsin 


,r h 


rindex t 


( 5 . 3 - 9 ) 


where RINDEXj^ is the refractive index in layer (L). The derivative dz/dy 
can be isolated and simplified by appropriate series expansions. It is used 
in the code to develop the true path as a function of the nominal path, and 
so to provide the correct boundaries, etc., at which to evaluate (TTS) and 
(TTD). 

Scattering Factor 

Besides the average attenuation discussed in the preceding subsection, 
the contribution to the output from a path increment (DD) depends on the 
previously introduced scattering factor of the form 

T, p ± (T|») DD 
i 

Since the atmosphere being modeled contains both Rayleigh and aerosol 
scatterers, there are two a ' s and p's that have to be evaluated. In both 
cases, the ( cr ) is the local extinction coefficient, denoted as (RDTAU) for 
Rayleigh scatterers, and (ADTAU) for aerosol scatterers. Also in both 
cases, the polarized output (if any exists) can be expected to be linearly 
polarized in the plane perpendicular to the scattering plane, so that in both 
cases two polarization parameters of the phase function constitute sufficient 
input data. Finally, for both cases it is essential to specify the scattering 
angle (</ 0 . Actually, this is done by specifying the cosine of that angle 
(COPSI). We simply form the dot product of a vector representing the 
input sunlight and a vector representing the light traveling to the detector. 

In the case of Rayleigh scattering, the phase function can be calculated 
analytically within the program. Ideally, the two polarization parameters 
are 6/1 67T and 6 COPSI / 1 6 77 . Sometimes it is appropriate to modify this 
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slightly. But for aerosol scattering, the situation is much more complicated. 
At present, we read in the phase function for selected values of (COPSI), 
usually 46 in number, and interpolate on the basis of (COPSI). 

Albedo 


Scattering from the atmosphere is only one of two possible mechanisms 
whereby light reaches the detector. There is also scattering from the 
surface of the earth, or some cloud layer above it. The total fraction of 
the incident light scattered from such a surface is its albedo. The simplest, 
and most commonly used angular dependence for surface scattering is the 
well-known Lambert law: 

y o 

radiance out = (irradiance in) ALBEDO — (5.3-10) 

where (/ j q ) is the cosine of the incidence angle. It is known that natural 
surfaces usually are not Lambertian. Nevertheless, better models are 
only occasionally available. Our model has incorporated a Feugelson law 
for single scattering from clouds, but otherwise uses Lambert's Law. 
The Feugelson law for the approximate albedo (A £ ) of a cloud of optical 
thickness (r) is given by 


A c “ 1 " 12 


2 ^o + 1 

1 + ( 1-^0 T 


7 ci y Q /3 + c 2 (3y Q - l)/8 {1 - 2y o (c x /3 - 1) } 
24 1 + (l-^-)T 


(5.3-11) 


where c^ = 2.48 and Cg = 3.70. 


5.3.3 Forward Scattering on the Sun Path 


After single scattering, the next less simple mechanism whereby light 
reaches the detector includes one or more nearly forward scatters on the 
path from the sun to the primary scatter. This subsection details the way 
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in which this contribution is handled by adapting the methods already 
described in the preceding sections. That is, we discuss modifications to 
the geometry, the average attenuation, the scattering fraction, and the albedo 
calculation. 

Geometry 

The geometry used to calculate the contribution to the output that 

includes some forward scattering on the sun path differs from the single 

scattering geometry only in that what is assumed to propagate along the 

sun path is a stream rather than a plane wave. The explicit definition for 

2 

the stream (Whitney, 1972) is radiance integrated against a cos weighting 
function. There occurs in the calculation a factor (FPRL), which is part 
of the stream definition, and factors (NORM) and (ESTPRL) that relate to 
the inversion of streams back to radiances at the output. This inversion 
is discussed in detail in a subsequent subsection on full multiple scattering. 

Attenuation 

The notion of a stream rather than a plane wave propagating along 
the sun path modifies the average attenuation in two ways. First, since 
forward scatters do not remove anything from the stream, forward scatters 
have to be excluded from the optical depths along the sun path. This is 
done by modifying the optical depths by a factor of (1-PHS), where (PHS) 
is a combination of streamlike phase function integrals that represents 
forward scattering from a stream into itself. The construction of (PHS) 
is detailed in the above mentioned paper. The modified optical depth from 
the sun is (MSTTS). With that replacement, the analog of (DAMP), the 
average attenuation, is 

MS DAMP = -Ae“ MSTTS " TTi yA(MSTTS + TTD) (5.3-12) 
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The second modification required is that we exclude from (MSDAMP) 
the case in which absolutely no scatters occur on the sun path, since we 
have already counted single scattering. This done by decreasing (MSDAMP) 
by that part of the average that represents just single scattering. Thus 
we put 


MSDAMP = MSDAMP - DAMP 


Scattering 

The scattering factor is also affected by the inclusion of forward 
scattering along the sun path. Where before it was appropriate to have a 
plane-wave transformation, or phase function p(0), it is now appropriate 
to have a stream transformation, or (PHS), evaluated for the nominal 
scattering angle (0). For both Rayleigh and aerosol scatterers, these PHS's 
are entered as data, precalculated at the angles that occur in the dodecahedral 
arrangement of streams introduced in (Whitney, 1972) namely 0°, 180°, 
63° and 117°. An interpolation routine then obtains the (PHS) for the desired 
angle (0). 

Albedo 

The final modifications required to accommodate a stream on the 
sun path concern the albedo calculations. For a Lambertian surface, a 
discretized version of Lambert's law is required. The continuous law 
involves the factor | cos(0)| In-, where 

tt = f | cos (0) | dft 
hemisphere 
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That is, the law is formulated to say that uniform radiance input provides 
uniform radiance output, just scaled down by the albedo. The discrete law 
must be formulated to say the same thing, but with the integral replaced 
by a sum over the dodecahedron angles, and the |cos(0)l replaced by 
averages of |cos(0)| with quadratic weighting functions. The normalization 
for all the averages is 

A A « 

(p*k) ^ da = 2 tt/3 (5.3-13) 

hemisphere 

The actual averages are as follows: for 0°, 

(z *k) 3 dfi = 3/4 (5.3-14) 

hemisphere 


< cos (6) 


.._a f 

2, ; 
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Thus the sum analogous to the integral is 

5 terms 

3/4 + ^ 9/20 = 3 (5.3-16) 


Thus for the discrete case, the (77) in Lambert's law is replaced by a 3. 
5.3.4 Forward Scattering on the Detector Path . 


We suppose now that the photons sustain one or more nearly forward 
scatters on their path from the primary scatter to the detector. This case 
is slightly more complicated than that ..discussed in the preceding section, 
and so entails a few more modifications, in regard to geometry, average 
attenuation, and scattering factor. 

Geometry 

The major geometric consideration with forward scattering on the 
detector path is the possibility of transfer of radiation from one line of 
sight to another at a slightly different angle. To handle this, we formulate 
for each line of sight an approximation for the probability that forward 
scattering ultimately contributes to that particular line of sight. This 
probability varies from zero for the uppermost line of sight to unity for a 
scan tangent to the earth. In between, the probability is approximated as 
follows: the difficult feature requiring approximation is that where forward 

scatters are concerned we must think not of contributions from a path 
increment, but rather from whole three dimensional regions. That is, if 
subsequent forward scatters are allowed, contributions to a given scan line 
come not from the line, but from a region surrounding the line. We formulate 
the above mentioned probability function as a line integral of the form 

(5.3-17) 
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where the path over which the integral is done traverses, roughly, the optical 

region that can, with no mechanism other than subsequent forward scattering, 

th 

contribute to the m scan line. The path proceeds from the detector, 

j. i_ 

through atmospheric layers down to and including the m L1 one, and then 
out again. The A r 1 to be associated with any particular layer has to be 
approximated in a rather gross way, because with the curvature of the 
earth, the different propagation directions in a stream can traverse 
drastically different optical depths. We begin by formulating Ar' in such 
a way that if the whole atmosphere were condensed in one uniform layer 
(whose geometric thickness would therefore be the atmospheric scale 
height), that layer would have the, optical depth tangent through itself: 

Ax' = EXTINCTION /2 RE SH (5.3-18) 

Thus for a layer of thickness KM, we put 

Ax' = EXTINCTION /2 RE KM /KM/SH (5.3-19) 

This formulation guarantees that for very small (KM), (r 1 ) converges to a 
well-defined function of the tangent height associated' with (PM) and is 
insensitive to the number and geometric thickness of layers that happen to 
be used in the atmospheric model. Finally, we must consider the effect of 
the scattering phase function. Clearly in the case of extreme forward 
scattering, the contributing region that Ar 1 represents should shrink to 
zero; that is, forward scattering becomes ineffective for transfer from 
one line of sight to another. We represent this volume effect by modifying 
Ar 1 to Ar' (1-PHS) 3 . 

Attenuation 

The attenuation for contributions to the output that will undergo one 
or more forward scatters on the detector path is modified in several ways. 
First, the attenuation on the detector path is modified by a factor (1-PHS) 
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to omit forward scattering just as the attenuation along the sun path was. 
The modified optical depth on the detector path is (MSTTD), so that at 
each point the attenuation is 

-MSTTS - MSTTD 
e 

Second, as we have already counted single scattering and forward scattering 

on the sun path, we have to subtract the portion of the function representing 
■MSTTD- TTD 

those events: e . Third, and novel to the case where there 

are to be subsequent forward scatters, there is no one unique and readily 
identified path over which this function should be evaluated and averaged. 
Recall that instead of contributions from a path increment, we are considering 
contributions from a spatial region. Instead of attempting to average the 
attenuation over a region, we simply leave it as a point function, evaluated 
at representative points (e.g. the tangent points of scan lines directed to 
the bottom of each atmospheric layer). 

Scattering 

The absence of a particular path to integrate over is again evident 
in the scattering factor, for no particular path length (DD) appears. Instead, 
(DD) is replaced by a standard length for subsequent forward scattering, 
which is just the /2 RE KM /KM/SH previously introduced in connection 

with the PM's. 

5.3.5 Full Multiple Scattering 

The final and most complicated mechanism whereby light reaches 
the horizon profile detector involves two or more scatters that cannot be 
classified as nearly forward. This section discusses the application of a 
method that handles these events, and also incidentally provides some 
additional information concerning the other mechanisms already discussed. 
The method is basically that introduced in Whitney, 1972. In that paper, 
the integro-differential equation of radiative transfer is discretized as a 
set of twelve coupled differential equations for streams. For symmetry 
the streams are arranged pointing to the faces of a regular dodecahedron. 
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Geometry 

The geometry of the horizon scan situation is modified slightly to 
achieve a speedy application of the method without significant loss of 
accuracy. Specifically, thicker layers are used (usually 20 km), and they 
are imagined to be flat until the final step, where curvature is introduced 
as a correction. Because the layers are flat, the introduction of a cloud 
is simplified; the analog of 


SQ 


/KM - /SOLIDTIP 
/KM 


(5.3-20) 


becomes simply KM - SOLIDTIP 

KM 


Also because the layers are flat, the original continuous equation of radiative 
transfer can be expressed with inverse cosines of incidence singles. In 
the discrete case, these are all replaced by inverses of cosines averaged 
with quadratic weighting functions; that is, inverses of the values 3/4 and 
9 /20 previously derived in connection with the discrete analog of Lambert's 
Law. 

Attenuation 

The attenuation of contributions to the output is calculated in much 
the same way as it was in the preceding section. All optical depths are 
scaled by a factor of the form (1-PHS), in order to eliminate forward 
scattering from the calculation. Then the integrals are done exactly. The 
general form of the integral required is 

lull 
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where r^, r^, 7 ^ . . . are the optical depths at which the first, second, third 
. , . scatters occur, and Cq, Cp c 2 ... are the cosines for the stream 
traveling from the sun, from the first scatter, from the second scatter, . . . 
. . . The various r's are constrained so that each of the exponentials is 
less than unity. The constraints are handled by a tree structure, or set of 
nested do-loops, that governs the calculations. The program first chooses 
the layer for the first interaction (Ll), then the inclination for the scattered 
stream (PI), then its direction (Ql), and then the layer for the second scatter 
(L2), and so on. If the first inclination (PI) is upward (downward) the 
second layer (L2) is required not to be below (above) the first layer (Ll). 

In the code, the multiple attenuation integral is not done all at once, 
but in pieces; since the 7-3 integral is done in a loop which is within the 
loop that does the integral, and which is in turn inside the loop that 
does the integral. As an example, the ^ integral is just 


Tid/Co - 1/ Ci) , 

e 

c i 

Because it has been isolated from the functions of r^, it is no longer true 
that the exponential is necessarily less than unity; it may as well be greater. 
However, when all the required integrals are multiplied together, the result 
is less than unity, as is physically required. 

A small novelty is required to handle the case where two or more 
scatters occur in sequence within the same layer. The problem for that 
layer is then a micro version of the problem for the atmosphere as a whole. 
We could handle it the same as the full atmosphere, breaking the layer 
into smaller sub layers. But such a procedure would be time consuming, 
so we attempt to get an approximate answer without resorting to it. For 
the case where all the scatters from (i) through (i+n) occur in one layer, we 
approximate the required integral with the product of the' T f T i + 1 * * * * T i+n 
integrals, but with an important modification. Consider the case where 
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all the inclinations are the same. The problem then reduces to a familiar 
one in Poisson statistics, and the required integral is known to be simply 
the above product, but divided by (n.'). Thus we induce that generally a 
factor of (l In!) should be incorporated. 


Scattering 


: We come now to the actual scattering operations for full multiple 
scattering. These are constructed by a procedure whose goal is to minimize 
the storage space and computational burden required. The main source of 
difficulty to be overcome is the inclusion of polarization, requiring four 
Stokes parameters Sq, s^, Sg. to characterize a stream and, in general, 
sixteen parameters to characterize a stream transformation. The input 
stream transformations are provided for streams in the Z-X plane separated 
by the dodecahedral scattering angles 0°, 63°, 117°, and 180°. For these 
cases, symmetry reduces the sixteen parameters to far fewer, so that storage 
space is minimized 


In particular, for 63° and 117°, there are only eight 




and another operating on 


and for 0° and 180° there are only three parameters, only two of 


parameters: one 2x2 matrix operating on 

CSfc 

which are independent. For any particular scatter that occurs in the code, 
the stream transformation is synthesized by combining the minimal stored 
iriformation with appropriate rotations of the plane-polarization parameters 
( s 2 ) • That is, first a 2x2 rotation is applied so that the polarization 
state in the nominal scattering plane acquires the name 1. Then the 
appropriate stream transformation is applied. Then another 2x2 rotation 
is applied to restore the polarization state in the nominal scattering plane 
to its final correct name. 


Program Control 

In the case of multiple scattering it is not enough simply to attenuate 
and scatter. The program must be prevented from expending excessive 
effort on negligible contributions to the output, and must be prevented from 
reconsidering the contributions it has more accurately calculated in the 
single scattering part of the code. These goals are accomplished by several 
tests in the code. 
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Of primary importance is the so-called twig pruning, or (ENDTEST) 
for multiple scattering. This test compares the current contribution to 
the output with a numerical standard. The standard is based on the output 
that was obtained from single scattering, on the total number of terms of 
the particular scattering order that are left to consider, and on an input 
parameter expressing the programmer's judgement of how accurate he 
wants to be. If the test is failed, it is judged unprofitable to further pursue 
that particular branch of the multiple scattering tree. Were it not for this 
test, the program would run for outlandishly long times without appreciably 
changing the output. 

Inversion of Streams 

The final step in creating the horizon profile is to take the streams 

Sa provided by the multiple scattering code, and invert them to estimates 
P A 

of physical radiance, S(k). This problem is the inverse of that involved in 

A 

forming the stream S a from the radiance S(k). Since forming S* is analogous 

P a Pa 

to performing an integral transform on S(k), it follows that estimating S(k) 

is analogous to inverting the transform. Because we do the original 

transform only for power n = 2, and not arbitrary (n), it is possible to 

invert only approximately. This section discusses what is thought to be 

an optimum procedure for approximating S(k). 


We begin by considering a simple case where S(k) is assumed to be 
equal to S(-k), and to have the form 


6 terms 


s(k) = 2 S( P> (f(P*£) 2 ‘ j) 


(5.3-21) 


with this form, it can readily be verified that 

5 terms 


S = — 
fa £ 3 


S (p) + ^ | S(p') 


(5.3-22) 
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Referring to Whitney, 19 72, we see that the inverse of this relationship 
would be 

S(f>) = § 5 . 

As in Whitney, 1972, we have , to separate forward-traveling from 
backward traveling radiation by extending these simple six term 
relationships to twelve terms. For the numbers used, the analogous 
inverse is 


d verms 


i s 

8 


a 


3 £' 


(5.3-23) 


s(+£) 




terms 

v s +6 . 


5 terms 

+ 40 2 S -£' 


(5.3-24) 

This is the basic relationship currently used to invert the array of output 
streams in the code. 
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Fig. 5.3-1 Horizon Scan Geometry 



Fig. 5.3-2 Effect of Refraction 
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